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This  grant  has  supported  research  in  the  following  areas: 

•  Optimal  Control  of  Distributed  Systems 

•  Adaptive  Identification 

•  Adaptive  Control 

•  Robust  Control 

The  Publications  section  of  this  report  lists  six  PhD  dissertations  and  sixteen  research 
papers  produced  by  research  supported  in  whole  or  in  part  by  this  grant.  Most  of  these  pub¬ 
lications  have  been  provided  to  AFOSR  previously.  The  Appendix  to  this  report  contains 
five  of  the  research  papers,  which  are  representative  of  the  project. 


Optimal  Control  of  Distributed  Systems 

The  main  focus  of  this  component  of  the  research  is  approximation  theory  and  numeri¬ 
cal  methods  for  design  of  finite  dimensional  compensators  for  optimal  control  of  systems 
represented  by  linear  partial  and  functional  differential  equations.  The  primary  class  of 
applications  is  large  flexible  space  structures.  Papers  dealing  mainly  with  approximation 
theory  and  numerical  methods  for  optimal  control  of  distributed  systems  are  [1, 2, 3, 4].  The 
paper  [5]  has  substantial  results  on  approximation  of  digital  input/output  models,  as  well 
as  results  on  infinite  dimensional  system  theory  that  are  useful  in  adaptive  identification 
and  control  of  distributed  systems.  The  paper  [6],  which  deals  exclusively  with  stability 
theory  for  a  class  of  wave  equations,  is  motivated  by  certain  stability  issues  that  arise  in 
approximation  theory  for  control  of  equations  governing  many  fleiuble  structures. 

Adaptive  Identification 

This  segment  of  the  research  has  dealt  with  fast  methods  for  real-time  adaptive  identification 
and  prediction  of  systems  with  unknown  parameters  and  unknown  order.  The  class  of  such 
methods  to  which  this  research  has  contributed  most  is  least-squares  lattice  filters.  The 
paper  [7]  demonstrated  the  application  of  a  lattice  filter  for  adaptive  identification  and 
prediction  of  an  experimental  flexible  structure.  The  paper  [8]  and  the  PhD  dissertation 
[9]  developed  a  new  lattice  filter,  called  an  unwindowed  lattice,  which  achieves  much  faster 
convergence  to  parameter  estimates  and  accurate  prediction  than  prewindowed  lattices. 


Adaptive  Control 

This  part  of  the  research  has  concentrated  on  adaptive  control  and  tracking  problems  for 
flexible  structures  and  manipulators  with  flexible  links  and  joints.  The  papers  [10,  11, 12, 
13, 14]  and  the  PhD  dissertations  [15, 16, 17]  have  developed  adaptive  contrd  methods  that 
appear  to  have  wide  application  to  flexible  space  structures  and  high-speed  manipulators. 
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Robust  Control 


The  goal  of  this  part  of  the  research  has  been  to  develop  nomerically  efficient  methods  for 
determining  robustness  margins  for  control  systems  with  uncertain  parameters  and  for  de¬ 
signing  plants  and  controllers  that  have  desired  robustness  mar^ns.  For  parameter  values 
within  these  robustness  margins,  the  systems  are  guaranteed  to  be  asymptotically  sta¬ 
ble.  The  papers  [18,  19]  presented  new  robustness  margins  along  with  numerically  efficient 
methods  for  computing  these  mar^ns,  and  the  paper  [20]  presented  a  method  for  numerical 
design  of  robust  control  systems.  These  papers  were  based  on  the  PhD  dissertations  [21, 22]. 
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Approximation  in  Control  of  Thermoelastic  Systems 

J.S.  Gibson  * 

Mechanical,  Aerospace  and  Nuclear  Engineering 
University  of  California,  Los  Angeles  90024 
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University  of  Southern  California,  Los  Angeles  90089 

G.  Tao 

Department  of  Electrical  Engineering 
Washington  State  University,  Pullmam,  WA  99164 

June  1, 1991 

Abstract 

This  paper  develops  an  abstract  framework  for  analysis  and  approximation  of  linear  ther* 
moelastic  control  systems,  and  Cor  design  of  finite-dimensional  compensators.  The  thermoelastic 
systems  in  this  paper  consist  of  abstract  wave  and  dlffnsion  equations  oonpled  in  a  skew  self- 
adjoint  fashion.  Linear  semigronp  theory  is  used  to  establish  that  the  abstract  thermoelastic 
models  are  well  posed  and  to  prove  convergence  of  generic  approximation  schemes.  Open-loop 
uniform  exponential  stability  Cor  a  subclass  of  thermodastic  systems  is  proved  via  a  Lyapunov 
function.  An  example  involving  the  design  of  an  optimal  LQG  compensator  Cot  a  thermoelastic 
rod  illustrates  the  application  of  the  abstract  theory.  Results  of  an  extenmve  numerical  study, 
including  a  comparison  of  the  dosed-loop  performance  of  different  compensator  deagns,  are 
presented  and  discussed. 


*Thii  author  was  supported  by  AFOSR  Grant  87-0373. 
(This  author  was  supported  by  AFOSR  Grant  87-0354. 


1  Introduction 

The  traasfer  of  energy  betifeen  ito  mechanical  form  and  heat  generally  has  been  ignored  as  a  source 
of  both  structural  damping  and  excitation  in  the  vast  literature  on  ccntrol  of  flexible  structures. 
Only  a  few  recent  p^era  have  considered  control  of  thermoelastic  structures  [1,  2,  3,  4,  5],  [6,  7,  8, 
9].  However,  the  thermally  induced  vibrations  that  hampered  the  recently  launched  Hubble  space 
telescope  have  highlighted  the  coupling  between  mechanical  vibration  and  heat  transfer  and  the  need 
to  model  and  control  thermoelastic  phenomena  in  flexible  structures. 

This  paper  has  two  main  objectives;  first,  to  develop  a  theoretical  framework  for  analysis  and 
approximation  in  the  design  of  feedback  control  systems  for  a  broad  class  of  linear  thermoelastic 
systems;  second,  to  illustrate  the  application  of  the  theory  by  presenting  the  most  interesting  re¬ 
sults  from  an  extensive  numerical  study  of  LQG  optimal  control  of  a  thermoelastic  rod.  Both  the 
theory  and  the  example  focus  on  numerical  methods  and  convergence  analysis  for  the  design  of 
finite-dimensional  compensators  baaed  on  finite-dimensional  approximationa  of  distributed  noodels 
of  thermoelastic  systems. 

By  a  thermoelastic  system,  we  mean  an  abstract  wave  equation  coupled  in  a  skew  self-adjc^t 
fashion  with  a  diffusion  equation.  While  some  of  the  theory  developed  here  pertains  specifically  to 
problems  in  which  the  generalized  wave  equation  is  second-order  (in  time),  much  of  the  theory  applies 
to  a  broader  class  of  problems,  mcluding,  for  example,  problems  in  which  a  Schrodinger  equation 
is  coupled  with  a  diffusion  equation.  In  this  paper,  we  are  particularly  interested  in  second-order 
generalized  wave  equations  because  they  are  conunon  in  flexible  structures,  but  the  results  here 
that  allow  a  more  general  class  of  wave  equations  are  intended  to  apply  also  to  problems  such 
as  thermal  blooming  in  lasers  [10].  Although  the  theoretical  framework  developed  in  this  paper 
handles  a  wide  variety  of  thermoelastic  systems,  it  is  not  clear  whether  our  hypotheses  hold  for  the 
thermo-viscoelastic  systems  with  memory  studied  by  Bums  et  al.  [1,  2,  3]. 

Our  philosophy  in  the  abstract  formulation  of  thermoelastic  control  systems  in  Section  2  and  in 
the  approximation  theory  in  Section  4  is  to  base  the  results  on  hypotheses  that  require  as  little  as 
possible  beyond  conditions  that  normally  hold  for  the  individual  wave  and  diffusion  equations.  This 
means  that,  in  analjning  a  particular  application,  most  of  the  work  is  done  on  the  uncoupled  wave 
and  diffusion  equations,  and  the  work  required  to  couple  the  systems  is  minimised.  For  example, 
in  verifying  the  hypotheses  for  Theorem  4.6,  which  concerns  convergence  of  approximations  to  the 
open-loop  thermoelastic  system,  once  the  convergence  conditions  for  independent  approximations  to 
the  uncoupled  wave  sod  diffusion  equations  are  verified,  no  further  work  is  necessary  to  guarantee 
convergence  cf  the  approximations  to  the  thermoelastic  system  when  the  straightforward  Galerkin 
scheme  that  we  assume  far  approximating  the  coupling  operator  is  used. 

The  approach  to  compensator  design  in  Sections  3  and  4.1  of  this  paper  is  to  ^proximate  an 
ideal  infinite-dimensional  LQG  compensator  with  a  sequence  of  finite-dimensional  compensators. 
However,  the  abstract  formulation  oi  thermoelastic  control  systems  in  Section  2,  the  approximation 
and  convergence  theory  in  Section  4.2,  and  the  result  in  Section  h  on  open-loop  uniform  expcmential 
stability  should  be  useful  in  any  method  for  analyns  and  design  of  controllers  for  thermoelastic 
systems. 

An  important  issue  in  both  convergence  of  the  approximating  compensates  and  performance 
of  the  closed-loop  systems  is  uniform  exponential  stability  of  the  open-loop  thermoelastic  system. 
While  several  authors  [6, 11,  4,  12, 13]  have  proved  strong  stability  for  various  linear  and  nonlinear 
thermoelastic  systems,  few  results  ha^  been  published  on  uniform  exponential  stability.  A  result 
in  [12]  on  integrability  of  the  energy,  when  applied  to  the  linear  case,  yidds  uniform  exponential 
stability  for  thermoelastic  rods  with  certain  sets  of  boundary  conditions.  Also,  a  recent  eigenvalue 
analysis  in  [6]  jrields  uniform  exponential  stability  for  linear  themooelastic  rods  with  the  same  sets 
of  boundary  conditions  to  which  the  result  in  [12]  applies.  The  proof  of  our  Theorem  5.1  uses  a 
Lyapunov  function  to  establish  uniform  exponential  stability  for  a  large  class  of  linear  thermoelastic 
systems,  but  does  not  improve  on  the  results  in  [6]  and  [12]  for  the  rod.  The  results  in  (6,  12]  and 
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our  Section  5  do  not  apply  to  the  set  of  boundary  conditions  for  which  uniform  ejq>onentiaJ  stability 
has  been  proved  recently  in  [14]. 

In  Section  6,  we  apply  the  theory  developed  in  Sections  2-5  to  design  finite  dimensional  compen¬ 
sators  for  a  thermoelastic  rod.  We  present  numerical  results  for  the  functional  control  and  estimator 
gains  that  represent  the  compensators  graphically.  We  also  compare  the  closed-Ioc^  eigenvalues  pro¬ 
duced  by  three  of  the  finite-dimensional  compensators  based  on  different  damping  models.  These 
eigenvalues  were  obtained  from  simulations  in  which  each  ctxnpensator  was  omnected  to  a  model  of 
the  rod  with  dimension  significantly  higher  than  the  dimension  of  the  compensator.  This  compar¬ 
ison  illustrates  the  importance  of  modelling  even  very  light  thermoelastic  damping,  or  possibly  an 
artificial  viscous  equivalent,  if  no  stronger  damping  mechanism  is  present. 
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2  Abstract  Thermoelastic  Systems 


Throughout  this  paper,  H  or  Hj  (J  =  0,1,2)  will  be  a  Hilbert  space  with  inner  product  (-,  •}  or 
(•,  •)j  and  corresponding  induced  norm  |  ■  j  or  |  ■  (/.  Also,  V  or  will  be  a  reflexive  Banach  space 
with  norm  ||  •  ||  or  ||  ■  ||j.  The  continuous  dual  of  V  will  be  denoted  by  V,  and 

(2.1) 

will  mean  that  V  is  embedded  densely  and  continuously  in  H,  which  implies  that  H  is  embedded 
densely  and  continuously  in  V  (see,  for  example,  [15,  16]).  In  this  case,  (-,  •)  will  denote  both  the 
ff -inner  product  and  the  duality  pairing  on  V  x  V. 


Lemma  2.1  Lti  V  and  H  he  related  as  in  (2.1),  lei  A  he  a  linear  isomorphism  (i.e.,  a  eoniinnons 
linear  hijection  with  continuous  inverse)  from  V  to  V  suck  that  A  is  dissipative  in  the  sense  that 

Re(v,Av)<0  \lveV,  (2.2) 


and  define 

Dom(A)  =  ^-»/f,  A  =  A\oom(A).  (2.3) 

Then  Dom(A)  is  dense  in  H  and  A~^  €  Also,  A  is  a  maximal  dissipative  operator  on  H. 

Proof  That  Oom(A)  is  dense  in  H  follows  &om  the  fact  that  H  is  dense  in  V  and  A~^  is  bounded 
from  V'  to  /f .  To  sm  that  A  is  maximal  dissipative,  suppose  that  there  exists  a  dissipative  linear 
operator  A  :  Dom(A)  C  H  H  that  is  a  proper  extension  of  A.  Since  7^(A)  =  H,  there  exists 
h  €  Dom(i4)\Dom(A)  and  v  €  Dom(A)  such  that  A  0,  Ah  =  0,  and  Av  =  h.  Then,  for  any  real  a, 
(w  -f-  ah,  A{v  -h  oh))  =  (v,  Av)  +  o[hp,  and,  for  sufficiently  large  o  >  0,  Re(t)  +  ah,  A(v  oh))  >  0, 
contradicting  the  dissipativity  of  A.  □ 


Theorem  2.2  Let  the  Hilbert  space  Hi,  the  reflexive  Banach  space  Vi  and  the  operator  Ai  he  as  in 
Lemma  2.1.  Let  the  Hilbert  space  and  the  reflexive  Banach  space  Vj  he  as  in  Lemma  2.1,  and 
let  At  he  a  linear  isomorphism  from  Vt  to  that  is  V^-eoemne;  i.e.,  there  exists  a  positive  real 
number  a  such  that 


Ba{,4,Ai4>)i  >  d  €  Vj. 

(2.4) 

Abo,  lei  £  €  B{Vi,Vi).  Define 

H  =  HixH2,  V  =  VixV3 

(2.5) 

and 

-[!■  :5;] 

(2.6) 

where  £'  €  B(V^a,  V^)  is  defined  hp 

rv)i  =  {*.£d>)2, 

(2.7) 

(i.e.,  C*  is  the  Bauach-apaee  adjoint  of  C).  Then  H,  V,  V  and  A  are  as  is  Lemma  2.1. 

Proof  Since  Ai  is  dissipative  and  At  is  Vj-coercive,  the  operator  (>4]  —  CAJ^C*)  E  B^Vt,  V^)  is 
Vj-coerdve.  Hence,  for  /i  €  V{  and  /]  €  the  pair  (viiVr)  €  V  given  by 


•>*(A,-£Ar*r)-*(£Ar‘/i-/»).  (2.8) 


is  the  unique  solution  to 


The  mapping  that  takes  (/i./j)  to  (vx.O])  is  clearly  bounded  from  V  s  V{  x  V,  to  V.  □ 


(2.9) 
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Remark  2.3  V/t  define  the  adjoint  operators  €  B(Vi,  Vj),  Aj  €  B(Vj,  V,)  and  A'  €  6(1^,  V') 
as  in  (2.7)  with  the  appropriate  duality  pairing  in  each  case.  Under  the  hypotheses  o/ Theorem  2.2, 
Al,  Al  and  A*  have  the  same  properities,  respectively,  asAi,  Aj  and  A. 

Remark  2.4  We  define  the  operator  L  :  Dom(I>)  C  Hi  Ha  to  be  the  restriction  of  C,  to  Dom(L)  = 
{^  €  Vi  :  Ctl>  €  Ha).  IfDom{L)  is  dense  in  Vi,  we  define  the  operator  L*  :  Dom(L*)  C.  Ha  —*  Hi 
to  be  the  Hilbert  space  adjoint  of  L  with  respect  to  the  Hi  and  Ha  inner  products.  It  can  be  shown 
that  I*  is  the  restriction  of  C*  to  Dom(£«*)  =  {^  €  Vj  :  €  Hi). 


For  the  class  of  systems  of  primary  interest  in  this  paper,  there  exist  Hilbert  spaces  Ho  and  Ha 
and  reflexive  Banach  spaces  Vo  and  Va  such  that  Vq  Ho  *-*  Vg  and  Va  •-*  Ha  *-*  V,  (with  each 
injection  continuous  and  dense).  The  thennoelastic  evolution  equations  have  the  form 

iv(t)  +  Pou»(0  +  >to«;(t)  +  £Stf(0  =  /o(0.  t  >  0,  (2.10) 

«(0  +  AaHt)  ~  Cowit)  =  fait),  t  >  0,  (2.11) 

where  Vo,Ao  €  B(Vo,  Vg').  £o  €  B(Vo,  Vj'),  Aa  €  B(Va,  V^),  €  £i(0,f;  B.)  for  t  =  0,2  and  all  f  >  0. 

We  assume  that  ^g  is  symmetric  in  the  sense  that 

i^,Ao4)o  =  {4>,Aot(>)Q,  ^,V’€Vo,  (2.12) 

and  that  Ao  is  Vo-coercive  and  Aa  is  Vj-coercive.  We  assume  that  27o  is  nonnegative  in  the  sense 
that 

Re(ti-,DoV')o  >  0,  V>  €Vo.  (2.13) 

To  derive  a  semigroup  generator  for  the  thermoelastic  system  in  (2.10)  and  (2.11),  we  first 
consider  the  semigroup  generator  corresponding  to  (2.10)  for  the  case  £o  =  0.  We  make  Vo  mto  a 
Hilbert  space  by  defining 

(tfr.^)v',  =  (V'Mo^)o.  ^.!^6Vo.  (2.14) 

Our  hypotheses  on  ^o  imply  that  the  norm  induced  by  the  inner  product  in  (2.14)  is  equivalent  to 
the  original  Vo  norm.  We  define 

Hi  =  Vox  Ho,  Vi  =  Vo  X  Vo,  (2.15) 

and  we  identify  Vo  with  Vq  in  the  first  c<xnponent  of  Hi  and  Vj  and  write  V/  =  Vo  x  Vg'.  It  follows 
that  Vi  /fi  ^  V{. 

Next  we  define 

^■=[-X  -p. 

That  .4i  is  an  tsomorphism  from  Vi  to  V{  follows  from 

■^'']€B(V;,V0.  (2.17) 

We  define  i4i  by  (2.3)  with  A,  A  and  H  replaced  by  v4i,  Ai  and  Hi,  respectively.  According  to 
Lemma  2.1,  Ai  generates  a  cootractioa  aeznig:oup  on  Hi.  (See  [15, 16, 17, 18]  for  similar  approaches 
to  obtaining  semigroup  generators  of  the  form  in  (2.16).)  Abo,  we  note  that  the  restriction  of  ~Aa 
to  Aa^Ha  generates  a  uniformly  exponentially  stable  analytic  contraction  semigroup  on  Ha-  For  the 
thermoelastic  system,  we  define 

r  =  l0  Co]£BiVi.Vi)  (2.18) 
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to  obtain  the  situation  in  Theorem  2.2  with  Ai  defined  by  (2.16).  The  corresponding  A  defined  by 
(2.6)  is 


A  = 


0/0 
-Ao  -Vo  -Ci 

0  £o  — Ai 


€  B(V,V') 


(2.19) 


where 


V  =  VoxVoxV3*-^H  =  VoxHoxHi^V'  =  VoxV^x  V,'.  (2.20) 

The  semigroup  generator  A  tor  the  thermoelastic  system  in  (2.10)  and  (2.11)  then  is  defined  by  (2.3). 
Explicitly,  the  domain  of  this  semigroup  generator  is 


Dom(A)  =  {(d.t&.fl)  €  V:A(^.M  €  H}.  (2.21) 

The  system  in  (2.10)  and  (2.11)  now  can  be  written  as 


x(0  =  Ai(0 +  /(<).  <>0,  (2.22) 

where  x(i)  =  ^  and  /  =  (0,/o,/j)  €  lj(0,f;  H)  for  all  f  >  0.  If  {r(<)  :  t  >  0} 

is  the  semigroup  generated  by  A,  the  mi/d  toluiion  to  the  initial  value  problem  consisting  of  (2.22) 
and  an  initial  condition  z(0)  =  (u;(O),ti)(O),0(O))  €  If  is 

z(t)  =  T(t)z(0)  +  f  T{t  -  s)/(s)ds,  t  >  0.  (2.23) 

Jo 
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3  The  LQG  Optimal  Control  Problem 

In  the  abstract  thenno«Iastic  system  (2.10)-{2.11),  we  consider  inputs  of  the  form 

=  +  t>0  (3.1) 

and  an  output  given  by 

y(t)  =  Cx(0  +  t^<).  <>0,  (3.2) 

where  z  is  the  mild  solution  to  (2.22),  u(<)  €  iZ",  7(0  €  iZ*,  y(0  €  IV,  u(t)  €  IP,  B  € 

B  €  and  C  €  B(H,IiF).  Also,  j  and  u  are  stationary  lero-mean  Gaussian  white  noise 

processes  with  covariance  matrices  T  and  R,  respectively,  and  R  is  positive  definite. 

The  linear-quadratic-Gaussian  (LQG)  optimal  control  problem  is;  given  the  output  y  in  (3.2), 
choose  u  to  minimize 

J(u)  =  Urn  f\{Qzit),zit))  +  u{ifRuit)]dt}  (3.3) 

I/-.00  tf  Jo 

where  Q  €  B{H,  H)  and  R  €  aie  self-adjoint  with  Q  nonnegative  and  R  positive  definite;  as 

in  (3.2),  z  is  the  mild  solution  to  the  thermoelastic  system  (2.10)-(2.11)  (or,  equivalently,  (2.22)) 
for  the  input  of  the  form  (3.1). 

In  view  of  (2.10)  and  (2.11),  the  operator  B  has  the  form 


where 

Bi  =  . . .  4,m]i  y  =  1, 2, . . . m,  {  =  0,2.  (3.5) 

The  operator  B  has  the  same  form.  The  operator  C  in  (3.2)  has  the  form 

C  —  [C3bi  Qa  Ci],  (3.6) 

where  Coi  €  B(Vo,  RF),  Cm  €  5(ffo.  /Z^)  and  C»  €  B{H,,  RF). 

Theory  for  the  infinite  dimensional  LQG  optimal  control  problem  with  bounded  input  and  output 
operators  can  be  found  in  [19,  20, 21, 22, 18].  We  triefiy  Bummarize  the  relevant  results  and  essential 
features  of  the  theory  here.  As  in  finite  dimensions,  the  LQG  problem  separates  into  a  deterministic 
linear-quadratic  regulator  problem  oo  the  infinite  interval  and  a  dual  state  estimator,  or  filtering, 
problem. 

First,  we  consider  the  regulator  problem,  which  is  to  choose  the  control  u  to  minimii^  the  integral 
in  (3.3)  when  both  noise  processes  in  (3.1)  and  (3.2)  are  sero,  the  output  operator  C  is  the  identity, 
and  tf  =oo.  If  thj  operator  pair  (A,  B)  is  unifonniy  exponentially  stabOizable  (i.e.,  there  exists  a 
bounded  linear  opttztor  K  such  that  A  -  BK  generates  a  uniformly  exponentially  stable  semigroup 
on  H)  and  the  pair  (Q,A)  is  uniformly  exponentially  detectable  (Le.,  the  pair  (A*,Q)  is  uniformly 
exponentially  stabilizable),  then  there  exists  a  unique  nonnegative  self-a4joint  solution  II  C  B(H,  H) 
to  the  operator  algebraic  Rkcati  equation 

A*n+nA-nBir*.B*n+Qso,  (3.7) 

with  II(Dom(A))  C  Dom(A*).  The  optimal  contrd  for  the  infinite-time  linear-quadratic  regulator 
problem  has  the  feedback  form 

o(f)  =  -Kzit),  t  >  0,  (3.8) 

where 

K  =  Rr^B*li^B{H,ir). 


(3.9) 


For  the  filtering  problem,  we  define 


Q  =  BTB’.  (3.10) 

If  the  pair  (C,  A)  is  uniformly  exponentially  detectable  and  the  pair  (A,  Q)  is  uniformly  exponentially 
stabilizable,  the  operator  algebraic  Riccati  equation 

+  EA*  -  tCTR-^CU  +  Q  =  0  (3.11) 

admits  a  unique  nonnegative  self-a4joint  solution  0  €  with  n(Dom(74*))  C  Dom(i4).  The 

minimum- variance  estimate  of  x(t)  given  v(r)  (r  <  t)  is  a  mild  solution  z(<)  to  the  evolution  equation 

i(0  =  AUt)  +  Bu(t)  +  k{v(t)  -  Ci(t))  (3.12) 

where 

k  =  ter  BT^e  BiR?,H).  (3.13) 

The  optimal  LQG  compensator  consists  of  the  filter,  or  state  estimator,  in  (3.12)  and  the  control 

Isiw 

u{t)  =  -Ki{t).  t>Q,  (3.14) 

with  the  control  and  filter  gain  operators  given  by  (3.9)  and  (3.13),  respectively. 

The  optimal  closed-loop  system  then  takes  the  form 

^(0  =  Sci(t  -  s)z(s),  0  <  s  <  t  (3.15) 

where  z(t)  =  (z(t),x(t))  ^  Z  =  H  x  H  and  {5ci(t)  :  <  >  0}  is  the  Co-semigroup  of  bounded  linear 
operators  on  Z  with  infinitesimal  generator 

^‘‘~{kC  A-  BK^  kC  ]  ’  ^  Dom(i4).  (3.16) 

If  {5(f) :  f  >  0}  and  {5(f) :  f  >  0}  are  the  seimgroups  df  bounded  linear  operators  generated  on 
H  by  infinitesimal  generators  A  —  BK  and  A  ->  kC,  respectively,  then  it  is  easy  to  show  that 

e(f)  =  5(f)e(0).  <>0,  (3.17) 

where  e(f)  =  z(f)  —  z(f).  Moreover,  if  for  some  real  a  and  M, 

||5(f)||  <  Afe— .  t  >  0,  (3.18) 

11^(011  <  i  >  0.  (3.19) 

then  for  each  b  <  a,  there  exists  a  omstant  Md  >  0  for  which 

l|5d(<)ll<A^.ie”‘*.<>0.  (3.20) 

Finally,  as  in  the  finite  dimensional  ease,  it  can  be  shown  that 

d(Ad)  =  <r(A  -  BK)  U  <r(v4  -  kC)  (3.21) 

where  c(Aei)  denotes  the  spectrum  of  the  closed-kx^  semigroup  generator  in  (3.16). 

We  note  that  the  uniform  exponential  stabilisability  and  detectability  conditions  stated  in  this 
section  are  sufficient  for  the  existence  of  unique  nonnegative  self-adjoint  solutions  to  the  operator 
algebraic  Riccati  equations  (3.7)  and  (3.11).  These  conditions  are  not  necessary  for  some  problems 
with  finite  rank  Q  and  Q.  A  sufficient  and  usually  necessary  condition  for  uniform  exponential 
stabilizability  and  detectability  is  that  the  open-loop  system  be  uniformly  exponentially  stable  except 
possibly  on  a  controllable  and  observable  finite-dimensional  subspaee. 
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It  ia  convenient  to  note  that,  since  C  iT"  and  Dom(A’)  =  RF,  there  exist  k  =  (ki, ....  k„) 
and  k  =  (ki,..  .,kf)  with  kj  and  kj  in  H  such  that 

[Kx]i  =  {x,kj),  x£H,  i=1.2 . m,  (3.22) 


kr  =  '^ kjTj  =^[kikj  ...  i^Jr,  r  €  /? .  (3.23) 

i=i 

Also,  €  H  implies  that  kj  =  (kj,i,kj,i,kjfl)  and  ij  =  (*i,i,ii,a,ij^)  with  €  Vo, 

€  i/oi  and  kj^ikj^  €  Hj.  It  follows  that 

{*•  ^i)  =  (^i  +  (t&t  4i.a)o  +  (^i  *>,3)3  (3-24) 

for  X  s:  e  H.  The  vectors  kj  and  kj  and  their  components,  kj^{  and  kj^i,  are  referred  to  as 

functional  control  and  estimator  (or  observer)  gains,  respectively. 
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4  Approximation  and  Convergence 

4.1  Approximation  Theory  for  the  LQG  Control  Problem 

An  approximation  and  convergence  theory  for  the  optimal  LQG  problem  for  iniinite-dimenaiona] 
systems  was  developed  in  [21,  23,  24,  18].  Here,  we  shall  first  briefly  summarize  the  generic  theory 
and  then  take  a  closer  look  at  it  in  the  context  of  abstract  thermoelastic  control  systems. 

Hypothesis  4.1  Then  txitU  a  se^eace  of  finHe-Hmenrionol  nbopaets  /T"  (n  =  1, 2, ...)  of  H,  sad 
scfseaces  of  openton  A"  €  Q"  €  B{H’',H''),Cr  € 

B{H’*,  BF).  The  openion  Q"  are  noMtgaiive  end  oelf-oijoini  for  etch  n. 

Ttom  here  on,  we  take  Q"  = 

Hypothesis  4.2  The  finite  dimeneional  etgehnie  Rieeati  efuations 

(A")*ip  +  IT  A"  -  +  Q"  =  o  (4.1) 

and 

A"n"  +  n^CA")*  -  n"(C")*.R-*c"n"  +  q"  =  0  (4.2) 

admit  unigue  nonnegative  telfiadjoint  solutions  II"  €  and  D"  €  B(H",H’*),  respeetivelp. 

We  define  gain  operators 

if"  =  i2-‘(B")*n"  €  Biir.ET),  (4J) 

and 

K”  =  n"(C")*B-‘  €  F(iP*,  IT),  (4.4) 

for  a  sequence  of  finite-dimensional  compensators  for  the  control  system  (2.10)-(2.11)  with  input  cd 
the  form  (3.1)  and  output  of  the  form  (3.2).  The  n-th  compensator  is  given  by 

«"(<)  =  -/ir"x"(<),  (4.5) 

i"(t)  =  A"i"(t)  +  B"«(<)  +  X"[y(t)  -  C"i"(t)l.  (4.6) 

The  resulting  closed-loop  system  is  then  given  by 

r"(<)  =  S3(<  -  s)z(s),  0  <  s  <  t  <  oo  (4.7) 

where  z"(t)  =  (i"(<),i"(t))  €  Z"  =  B  x  B",  and  {S]j(t) :  t  >  0}  is  the  Ce-semigroup  of  bounded 
linear  operators  on  Z"  with  infinitesimal  generator  A^  :  Dom(A^)  C  Z"  Z"  given  by 


=  [^ 


A  -BB" 

"C  (A"  -  B"B"  -  B7*C7"] 


]• 


Dom(AS)  =  Dom(A)  x  B". 


Since  B"  €  B(B",ir)  and  B"  €  B(Br,B")  we  have 

[B"x"l^  =  (tJ',x")j,,  y  =  l,2,”-,m 
for  x"  €  B"  and  ^ 

for  r€  with  i^,k^  €  B",  i  =  1,2,. ..,m,  j  ss  l,2,--*,p. 

The  convergence  theory  can  be  sunomarised  as  foOows.  We  will  refer  to  the  following  finite 
dimensional  semigroups: 


(4.8) 

(4.9) 
(4.10) 


T"(<)  =  e^‘’,  S"(<)  =  ^(1)  = 

and  their  adjoints  T"(t)*,  S"(t)*  and  ^(t)*. 


(4.11) 
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Hypothesis  4.3  For  taeh  n,  there  exists  s  linear  mapping  P"  from  H  onto  H"  sncA  that 


lim  P^z  —  *, 

x£H. 

(4.12) 

n>^oo 

For  each  z  £  H  and  each  <  >  0, 

hra  T'{t)P^z  = 

T(t)z, 

(4.13) 

lim  T'ltyP^z  = 

T(tyz. 

(4.14) 

n^oo 

where,  in  each  ease,  the  convergence  is  uniform  tn  t  fort  m  hounded  intervals.  Also, 

lim  =  Bu,  ue.BT', 

(4.15) 

n^co 

lim  =  Qz, 

z€H, 

(4.16) 

fl<^00 

and 

lim  CP^z  =  Cx, 

z€H. 

(4.17) 

n—oo 

If 

sup  ||n''||  <  oo  and 

fl 

8up||n"||  <  oo 

91 

(4.18) 

and  there  exist  positive  constants  M  and  a,  independent  of  n,  for  which 

|15"(1)||  <  Afe-*‘,  and  ||^(01l  <  <  >  0, 

(4.19) 

then  the  algebraic  Riccati  equations  (3.7)  and  (3.11)  admit  bounded  nonnegative  self-a4jomt  solu¬ 
tions  H  and  jQ,  and 


lim  n"P“x  =  n«, 

fl-»00 

*€/f. 

(4.20) 

lim  n"p'’i  =  n*, 

n«-*oo 

x£H. 

(4.21) 

Also, 

lim  S"(t)P"x  =  S(t)x, 

n-woo 

x£H, 

(4.22) 

and 

lim  ^{t)P”z  =  S{t)x, 

ti^OO 

z€H, 

(4.23) 

with  the  convergence  uniform  in  t  in  bounded  <-intervals.  If,  in  addition,  the  operators  Q”  and 
are  coercive  and  bounded  away  from  0  uniformly  in  n,  then  the  uniform  boundedness  of  ||II"||  and 
||n''||  yields  the  existence  of  positive  constants  M  and  a  independent  of  n  for  which  (4.19)  holds. 

The  easiest  way  to  guarantee  (4.18)  and  (4.19)  is  to  show  that  there  exist  positive  constants  Af 
and  a,  independent  of  n,  for  which 

|ir*(0||<Me— .  <>0,  (4.24) 

although  such  a  uniform  decay  rate  for  the  ^proximating  open>lo<^  semigroups  does  not  always 
exist.  When  (4.18)  holds  but  the  semigroups  {^(t)  :  I  >  0}  and  {^(f) :  (  >  0}  are  not  necessarily 
uniformly  exponentially  stable,  uniformly  in  n,  then  bounded  nonnegaiive  self-a4)oint  scdutions  II 
ud  n  to  (3.7)  and  (3.11)  exist,  but  1I«  and  RU  are  guaranteed  only  to  converge  weakly  to  II  and 
n,  respectively,  as  n  — » oo. 

When  the  strong  convergence  in  (4.20)  and  (4.21)  holds,  we  obtain 

^hm  Hif-P"  -  A'lliKH.it-)  =  0. 
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lim  ||K'"  —  ^||a(itr,i7)  =  0, 

#1^*00 

(4.26) 

and  therefore 

=  i  =  l,2,--,m, 

n-^oo  • 

(4.27) 

and 

Umkj=kj,  i  =  l,2,-”,p, 

fl^OO  • 

(4.28) 

in  17.  If  we  define  :  Z  —*  Z’*  hy 

i5=[J  ^]. 

(4.29) 

then  we  obtain  further  that 

53(<)P3z  =  S«j(0^,  i^Z, 

fl^OO 

(4.30) 

uniformly  on  bounded  t— intervals. 
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4.2  Abstract  Approximation  Theory  for  Linear  Thermoelastic  Systems 

Now  we  condiser  the  conatniction  of  the  approximating  finite  dimensional  subepaces  H’',  the  map¬ 
pings  P"  and  the  operators  A",  B".  Q",  etc.  We  establish  a  generic  approximation  theory  for 
abstract  linear  thermoelastic  systems  that  includes  relatively  easily  verified  sufficient  conditions  for 
the  convergence  in  Hypothesis  4.3. 

We  assume  the  hypotheses  of  Theorem  2.2. 

Hypothesis  4.4  For  j  s  1,2,  and  n  =  1,2,3,...,  H7  ts  o  finite  dimensional  anbtpace  ofVj  and 
Aj  €  anch  that  the  followng  eonditiona  hold. 

(i)  For  each  wj  €  0  =  1|2),  <Aere  exitU  a  sequence  o"  €  H”  suck  that 


(4.31) 

(ti)  For  each  n,  j4"  is  dissipative;  i.e., 

<  0.  v6H^. 

(4.32) 

(Hi)  For  each  f  €  V(  and  each  real  A  >  0, 

(A->lJ)->Pr7-^(A-^i)-V 

(4.33) 

and 

(A->in-‘pr/-i^(A->ti)-v, 

where  pp'  €  Hf)  is  defined  by 

(4.34) 

=  {v,f)i,  V  €  Hp,  j  =  1,2. 

(4.35) 

(tv)  7%ere  exsits  a  positive  constant  a  such  that,  for  alt  n. 

Re(v,  A^v}3  >  o||v||5,  V  €  Ba . 

(4.36) 

(r)  For  each  f  6  and  each  real  A  >  0, 

(A  +  A5)-*^/-S^(A  +  ^a)-V 

(4.37) 

and 

(A + )-‘^/  (A +>4;)-7- 

(428) 

Ramark  4.5  The  operator  Pp  restricts  s  functional  f  €  VJ  to  Hp  and  identifies  /|h-  with  an 
element  of  Hp  wia  the  Riesx  map  for  Bp.  Iff  esa  he  identified  with  an  element  of  Hi  (via  the  Eiesx 
map  for  Hj),  then  Ppf  is  the  Hi-  projeetion  of  f  onto  Hp. 

With  Pp  defined  by  (4.35),  we  define  I"  €  BiH^,HJ)  and  I."*  €  B{HS,H^)  by 

L’'s:Ppc\ai>  or  (os.X'"vi)3  =  {ps,£«i)s.  vi  €  BTiVi  €  B? ,  (4J9) 

L”’=Pp£*\Bf  or  (»i,£"*W3)i  =  wi€Br,t»j€BJ.  (4.40) 
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Hence  L"*  is  the  Hilbert-space  s4joint  of  X".  The  operator  I."  is  a  straightforward  Galerkin  ap¬ 
proximation  of  C.  On  the  other  hand,  Hypothesis  4.4  does  not  require  that  A"  and  be  Galerkin 
approximations.  Next,  we  define 


if"  =  Ifp  X  if  J 

(4.41) 

and 

(4.42) 

Theorem  4.6  For  fi  €  V{,  /]  €  V]  and  A  >  0, 

(4.43) 

and 

)  Ji-cA-xT'  ( 2 ) 

(4.44) 

Proof 

For  /i  €  V(  and  /j  €  V^,  we  set 

(4.45) 

and 

(4.46) 

We  note  that  (4.46)  is  equivalent  to 

t?  =  (A  -  .4j)-‘(pr/i  -  =  (A  -  A?)-*^(/,  -  rvj) 

(4.47) 

and 

uj  =  (A  -I-  Air\P^h  +  i-vj)  =  (A  +  ^5)"»P,"(/j  +  £vj). 

(4.48) 

Substituting  (4.47)  into  (4.48)  yields 

[(A  +  A^)  +  £"(A  -  .47)-‘l"*K  =  Pxh  +  4i"(A  -  X?)-*#!*/)- 

(4.49) 

ikom  (4.32),  (4.36),  (4.39)  and  (4.49),  we  have 

<»lh"ll’  <  Re((vM(A  +  >»5)  +  i"(A-Xr)-»I">?),) 

=  Re((«J,.^/,),  +  (t^,I-(A  -  A^r^Prfih) 

(44)0) 

=  Ile((vJ,/,),  +  (eJ.£(A  -  >4?)-*Pr/i>») 

<  ll»Jllj(l/>lv;  +  II^II  •  ll(A  - 

Since  (4.33)  implies  that  ||(A  —  i47)~*^/i||i  is  bounded  in  n,  (4.50)  shows  that  ||v]  Hi  is  bounded 
in  n.  Ilten,  it  follows  from  (4.33)  and  (4.47)  that  ||of  |||  is  bounded  in  n. 

Next,  we  note  that,  for  x  =  (xt,  zj)  €  IP', 

Re(z,  (A  -  A")z>  s  Re({z,.  (A  -  ilj)z,),  -I-  <z,,  (A  +  ^;)z,),)  >  a||z,|li  +  A|z|».  (4.51) 


We  set 


V 5?  ;■  V  (A+vt5)-»^(/,+£vo 


) 


(44i2) 
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and 


(4.53) 


^  =  ('1 

Then,  recalling  (4.47)  and  (4.52)  yields 

(z?,  (A  -  ^DzDx  =  (z?. (A  -  A7K  -  (A  - 

=  -(zr.r(«j  -  »,))!  =  -{zr, rzj)x  -  (zj,r(t;j  - «,)),  (4.54) 

and  similarly  (4.48)  gives 

(zj.  (A  +  >i5)zj),  =  (zJ.(A  +  A^)vq  -  (A  +  A;)t!5)j 

=  {zj.r(t.r  -  «x)),  =  (zj.£z?)a  +  (zJ.AC?  -  «»)),.  (4.55) 

Hence, 

Re((zr,(A  -  ADzDi  +  (z?.(A  +  >4;)z;),)  =  Re(-(zr.r(c;  -  Vi))i  +  (zJ.AC?  -  t;i)),).  (4.56) 
In  view  of  (4.51)  then, 

-  «Jllj  <  M\h  •  ll^ll  •  lie?  -  V7h  +  ll*?lb  •  ll^ll  •  lie?  -  villi  (4.57) 

According  to  (2.6),  (4.45),  (4.52)  and  conditions  (tit)  and  (v)  of  Hypothesis  4.4, 

lim  115"  —  willi  =  0  and  lim  ||5J  —  wjHj  s=  0.  (4.58) 

fl^OO  fl^OO 

Hence,  ||5"||  is  bounded  in  n,  and  we  have  seen  that  ||v"||  is  bounded  in  n.  Hence  ||z"||  is  bounded 
in  n.  Therefore,  (4.57)  and  (4.58)  show  that  O]  converges  in  Vj  to  «].  Then  (4.47),  condition  (nt) 
of  Hypothesis  4.4  and  (4.58)  show  that  v"  converges  in  Vi  to  vi,  and  (4.43)  ii  proven.  The  procf  of 
(4.44)  is  the  same  except  that  all  operators  except  and  ate  replaced  by  their  a<(joints.  □ 

Hypothesis  4.4  holds  for  most  common  approximation  schemes,  Gakrkin  schemes,  in  particular. 
The  foUowmg  theorem  establishes  conditions  (iv)  and  («)  of  Hypothesis  4.4  when  AJ  represents  a 
Galerkin  approximation  of  Aj. 

Theorem  4.7  Atattme  the  kfpotheaes  of  Thewem  2.2  rt§arding  Hi,  Vi  ond  Ai,  sad  sssame  con* 
dition  (i)  o/ Hypothesis  4.4  /or  j  =  2.  Define  A]  €  BiHgtEJ)  if 

AJ  =  P^Ai  In;  or  (v.A^to))  =  {v,Aiw)i,  »,to  €  HJ.  (4.59) 

Then  conditions  (ivj  and  (v)  o/ Hypothesis  4.4  hoH. 

Proof  condition  (iv)  is  immediate.  Tb  prove  (4J7),  let  /  €  V^'  and  set 

•  =  (A  +  v4,)-»/.  (4.60) 

••  =  (A  +  A;)-»^/.  (4.61) 

Also,  let  v"  €  if]  such  that  ||v"  —  «||3  converges  to  0.  Then 

o||«-|g  <  I(v-.(A  + A;)v"),|  =  \{v\PSf)i\  =  |(v-./),|  <  ||v-||,I/k..  (4.62) 

Hence  ||v"||2  is  bounded  in  n.  Next, 

o|tv-  -  5*11’  <  |(v-  -  5*.  (A  +  A;)(v-  -  5-)),| 

=  |(v"  -  5*. PJ/),  -  (v*  -  v-,(A  +  a;)5*),|  =  Kv*  -  5*. / -  (A  +  Aa)5-),1.  (4.63) 
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Since  v"  converges  in  V3  to  o  and  A)  €  B(V2,  V^),  it  follows  that  ||v''||3  is  bounded  in  n  and  (A-f><3)v" 
converges  in  Vj  to  (A-|><43)v  =  /.  Therefore,  (4.63)  shows  that  ijv”  —  v’'||3  converges  to  0  as  n  — » 00, 
so  that  o"  converges  in  Vs  to  «. 

The  proof  is  the  same  when  Aa  and  ^43  are  replaced  by  their  adjoints.  □ 

When  has  the  form  (2.16)  and  is  a  Gakrkin  approodmation  of  Ai,  condition  (tti)  of 
Hypothesis  4.4  can  be  proved  either  by  arguments  siinilar  to  the  proof  of  Theorem  4.7  or  by  projec¬ 
tion  arguments  like  those  in  [18].  Also,  see  [17]. 

Usually,  the  operator  P"  in  Hypothesis  4.3  is  the  P-projection  onto  P”  =  P"  x  *0  that 
condition  (i)  of  Hypothesis  4.4  guarantees  (4.12).  In  this  case,  if  €  Hj,  then  P^(fi,fa)  = 
(Pf  fi^Pafi)  (recall  Remark  4.5).  Hence,  it  follows  from  Theorem  4.6  and  the  Ikotter-Kato  theorem 
[25]  that  the  approximating  open-loop  semigroups  Tnif)  and  T^{t)  converge  as  in  Hypothesis  4J. 

Also,  when  P"  is  the  P^-projection,  it  is  naost  common  to  define  the  approximating  input, 
state-weighting  and  output  operators  by 


P"  =  P"P 

(4.64) 

g"  =  P"<3U.. 

(4.65) 

cr  =  c\H; 

(4.66) 

so  that  (4.15),  (4.16)  and  (4.17)  follow  from  (4.12). 


4.3  Matrix  Representations  of  Approximating  Operators 

We  assume  now  that  A"  has  the  form  in  (2.16),  that  C  has  the  form  in  (2.18)  and  that  Pi  and  Vi 
have  the  forms  in  (2.15).  Then  Pf  has  the  form  Hq  x  Pq  with  P^  C  V^.  We  assume  that,  for  each 
n,  Pq  is  the  span  of  a  finite  nun^er  of  basis  vectors  and  P3  is  the  span  of  a  finite  number  of 
basis  vectors  ef  j.  (The  spaces  Pg  and  P,  may  have  different  dimensions.) 

Also,  we  use  Galerkin  ^proximations  of  both  Ai  and  Aa-  The  matrix  representation  of  the 
operator  A"  in  (4.42)  is  then 


matrix  representation  of  A"  =  [A"]  = 


0  7  0 

Mn~^  trn  1/"”*  IT"  

0  /V0  ••AI0  — Al0  A3 

0 


(4.67) 


where 


=  K<.«Sj)o]  Af?  =  Ke5.r.e5j>»] 

=  (4.68) 

^  =  lKoT>o^M 

The  matrix  representation  of  the  operator  P”  in  (3.1)  and  (3.4)  is 


0 

[P-]=  Afo-'*K«S^.^)o]  . 

.  Afy*Ke5.i.*v>J  . 


(4.60) 


and  the  matrix  representation  of  the  operator  B"  is  sinular.  The  matrix  representation  of  the 
operator  C  in  (3.2)  and  (3.6)  is 


[C"]  =  [[Ci,ie;j(Co,eS.J[C,e;j]. 


(4.70) 


To  diacuos  the  matrix  representations  of  the  operators  Q”,  (f*,  n"  and  D",  it  k  convenient  to 
define  basis  vectors 


*0,»  ~  (*O,j>0|0)  *1,4  —  *2,4  “  (0>  Ot  *a,,-) 

and  the  block-diagonal  matrix 

A/"  =  diag{AfJ./fJ,AfJ}. 

The  matrix  representations  of  Q"  and  are 

m  =  (Q"]  =  A/"*‘Ke?..,.<5e7,j)].  i'./  =  0. 1,2. 


(4.71) 

(4.72) 

(4.73) 


The  matrix  representations  [11"]  and  [0"]  of  II"  and  H",  respectively,  are  determined  by  solving 
Riccati  matrix  equations  equivalent  to  the  operator  equations  (4.1)  and  (4.2).  The  form  [11"]  is 
like  that  of  [Q"],  and  in  general  neither  of  these  matrices  k  symmetric.  Hence,  rather  than  solving 
the  matrix  representation  of  (4.1)  directly,  it  k  preferable  to  premultiply  the  matrix  representation 
of  (4.1)  by  Af"  to  obtain  a  Riccati  matrix  equation  that  can  be  solved  for  the  symmetric  noatrix 
M"[n"].  Also,  instead  of  solving  the  matrix  representation  of  (4.2),  it  k  preferable  to  postmultiply 
the  matrix  representation  of  (4.2)  by  Af"  to  obtain  a  Riccati  matrix  equation  that  fjm  be  solved 
for  the  symmetric  matrix  [n"]A/"  (See  [18].) 

Finally,  it  follows  from  (4.3)  and  (4.4)  that  the  approximating  fimctional  ccmtrcd  and  estimator 
gains  in  (4.9)  and  (4.10)  are  given  by 


[kj  *;...*"]  =  e"A/"'‘[n"]Af"[R"]7r*,  (4.74) 

^  [k^ij  ...*;]  =  e"[n"]A/"'*[c"]’'R-S  (4.75) 

where 

*■"  =  I  (*7.4]  (*7.4]  1*7,41]  (4.76) 

and  [cq,,-],  for  example,  k  the  row  matrix  containing  the  bask  vectors  Cq,;  in  order.  See  [18]  for 
details  on  computing  similar  functional  gains. 
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5  Stability  of  the  Open-loop  System 

We  consider  the  system  in  (2.10)  sad  (2.11),  and  we  define 

Dom(i4o)  =  Ao  =  ^loomCita)-  (5-1) 

Since  >>  symmetric  and  Vo-coercive,  Ao  is  self-adjoint  and  Vr-coercive.  We  recall  the  operators 
£  and  £o  io  (2-18)  and  note  that  £o  €  B{Vo,  V^). 

In  this  section,  we  assume  that 

£o  =  Lo€B(Vo,Ho).  (5.2) 

and  we  assume  that  there  exists  a  positive  real  number  a  such  that 

I>om(i4o)  =  {o  €  Dom(£o) :  Lov  €  Dom(Lo)}  uid  Ao  =  oLlLo,  (5.3) 

where  Lq  is  the  Hilbert-space  a4joint  of  Lo  with  respect  to  the  Ho  •nd  Hi  inner  products  (recall 
Remark  2.4).  In  this  case, 

{v,w)vt=-a{Lov,Lov>)o,  t>,w€Vo.  (5.4) 

The  conditions  (5.2)  and  (5.3)  are  common  in  thermoelastic  structures  because  the  thermal  stress 
enters  the  equation  governing  mechanical  vibrations  in  the  same  way  as  the  stress  due  to  elastic 
deformation  [26,  27]. 


Theorem  5.1  Assume  the  conditions  stated  so  far  in  this  section  and  that  the  damping  operator 
X>o  is  symmetric  (in  the  sense  of  (2.12)).  If  the  range  of  the  operator 

s 

Ao  =  LoA^^  (5.5) 

is  in  Vi  or  if  Vo  is  Ho-coercive,  then  the  semigroup  generated  on  the  space  H  ta  (2.20)  ig  the 
operator  A  defined  in  (2.19)-(2.21)  is  uniformly  exponentially  stable. 


Proof  First  consider  the  ease  where  7S(Ao)  C  Vj  but  Po  i*  not  necessarily  ifo-coercive.  It  is 
clear  that  Ao  €  B(Ho,Hi)  and  a;  €  B(Hi,Ho).  Hence.  7^(Ao)  C  Vj  implies  Ao  €  B{Ho,Vi). 
Furthermore,  it  can  be  shown  that  AS  €  B(Hi,Vo)  and  aLohi  **  i^rprojection  onto  Ti{Lo). 
Now  define  the  following  self-adjoint  bounded  linear  operator  on  H: 


Q  = 


cl  0 

/  <rl  -2oA5 
0  — 2aAo  cl 


(5.6) 


where  o’  is  a  positive  real  number.  For  c  sufficiently  large,  Q  is  J7-ooerdve.  Also,  rince  TU^ho)  C 
and  7^(AS)  C  Vo,  QV  C  V.  For  Po  =  0  and  x  =  (o,fi,«)  €  Dom(A)  C  V, 

Re((?x,Ax)  =  Re<gx,>4x>  = 

(5.7) 

-I|f|l2  -  |AI2  -  cBjt{»,Ait)i  -I-  (20  -  l)Re(tf,Xor),  +  2o(LoA;tf,0*  +  2aIle(Ao.'i,/39),. 

Since 

|(Aofi,^30)]|  <  IIAofills  •  MxHbci'j.v')  •  ll^lb.  (5A) 

and  A«  €  B(Ho,  Vj),  it  follows  from  (5.7)  that,  for  c  sufficiently  Urge,  there  exists  a  positive  real 
number  fi  such  that 

Re(gx,  Ax)  <  — ^|x|’,  X  €  Dom(A).  (5.9) 
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When  T>o  ^  0,  the  right  side  of  (5.7)  has  more  terms,  but  (5.9)  can  be  obtained  in  a  siinilar  manner. 
The  generalized  Schwarz  inequality  |(v,I^ob)op  <  |(v,7>ot)o|  -  |(A,I)ob)o|  is  useful. 

If  2>o  »  /To-coercive,  then  replacing  a  with  0  in  (5.6)  allows  (5.9)  to  be  obtained  for  a  sufficiently 
large  and  some  positive  □ 

Remark  5.2  TAe  condition  72(Ae)  C  Vj  is  equivalent  to  the  following  two  conditions  combined: 

Dom(i;)  C  Vi  (5.10) 

and  tAere  exists  a  real  number  fi  such  that 

Ibllz  <  V  €  DoniLDn (5.11) 

Remark  5.3  To  generalise  Theorem  5.1  to  the  case  where  Dq  is  not  symmetric,  we  would  have  to 
impose  further  conditions  on  Do,  which  would  take  ns  begond  the  focus  of  this  paper. 


The  hypotheses  of  Theorem  5.1  hold  for  many  but  not  all  linear  thermoelastic  systems  that 
seem  likely  to  be  uniformly  exponentially  stable.  In  most  applications,  the  conditions  (5.10)  and 
(5.11)  restrict  the  combinations  of  boundary  conditions.  For  example,  if  (2.10)  and  (2.11)  represent  a 
thermoelastic  rod,  as  in  the  example  in  the  next  section,  (5.10)  and  (5.11)  hold  for  Dirichlet  boundary 
conditions  on  the  wave  equation  at  both  ends  of  the  rod  and  Nuemann  boundary  conditions  on  the 
heat  equation  at  both  ends,  and  for  various  other  combinations.  However,  (5.10)  and  (5.11)  do  not 
hold  for  Dirichlet  boundary  conditions  on  both  equations  at  both  ends  of  the  rod. 

Recently,  J.U.  Kim  [14]  has  proved  that  the  linear  thermoelastic  tod  with  all  Dirichlet  boundary 
conditicms  is  uniformly  exponential!  stable.  We  have  tried  without  success  to  modify  the  hypotheses 
of  Theorem  5.1  to  cover  this  case.  Also,  it  might  be  possible  to  apply  the  methods  used  by  Slemiod 
in  [12]  for  nonlinear  thermoelasticy  to  prove  uniform  exponential  sttd>ility  for  the  linear  all-Dirichlet 
case,  but  (5.10)  and  (5.11)  hold  for  the  boundary  conditions  treated  in  [12].  Hansen  [6]  has  shown 
that  all  of  the  eigenvalues  are  bounded  strictly  to  the  left  of  the  imaginary  axis  for  the  linear 
themooelastic  rod  with  all  Dirichlet  boundary  conditions,  but  Hansen’s  analysis  suggests  that  the 
eigenvectors  do  not  form  a  Riesz  basis. 

The  conditions  (5.10)  and  (5.11)  say  that  the  operator  As  in  the  diffusion  equation  is  bounded 
in  a  certain  sense  with  respect  to  the  stiffness  operator  Aq.  We  believe  that  some  such  relative 
boundedness  is  necessary  for  uniform  exponentially  stability.  A  numerical  ejq;>eriment  in  which 
we  used  the  one-dimensiunal  wave  equation  for  (2.10)  and  a  fourth-order  one-dimensional  partial 
differential  operator  for  As  in  (2.11)  yielded  a  sequence  of  complex  eigenvalues  that  appeared  to 
approach  the  imaginary  axis  asymptotically. 
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6  An  Example  and  Numerical  Results 

6.1  Linear  Model  of  a  Thermoelastic  Rod 

We  consider  the  axial  vibrations  of  a  visco-thermoelastic  rod  that  is  clamped  and  insulated  at  both 
ends.  The  length  of  the  tod  is  normalized  to  1.  Control  actuation  is  produced  by  a  single  force 
directed  parallel  to  the  rod  and  distributed  uniformly  over  the  rod  segment  171  <  >7  <  ih-  A.  sensor 
measures  axial  displacement  at  17  =  171  (i.e.,  the  left  end  of  the  rod  segment  over  which  the  actuator 
force  is  distributed).  Finally  we  assume  that  both  the  actuator  input  and  sensor  output  are  corrupted 
by  sero-mean  Gaussian  white  noise  with  unit  intensities. 

The  dynamics  of  the  plant  are  described  by  the  equations  at  one-dimensional  linear  thermoe- 
lasticity  (see,  for  example,  [26,  28,  13]),  which  consist  of  coupled  one  dimensional  wave  and  heat 
equations.  If  the  rod  has  Kelvin- Voigt  viscoelastic  damping  in  addition  to  thermoelastic  damping, 
then  the  state  equations,  boundary  conditims,  and  output  equation  are 


P^j(<.»7)  op(A  +  2/i)^^,^(t,i7)  + 

(6.1) 

80 

-I-Q'I,(3A  +  2p)—{t,  17)  =  6o(t7)u{t)  boivhit),  0  <  17  <  1, 

1  >  0, 

An  A^n 

PC^(<.»7)-«^(<,»7)  +  ^aL(3A-»-2/i)^^(<,i7)  =  0  0  <  17  <  1, 

t  >  0, 

(6.2) 

u;(t,0)  =  0  =  u>(<,l),  t>0. 

(6.3) 

|^(t,0)  =  0  =  |(t,l),  t>0. 

(6.4) 

y(t)  =  tu(<,f7i)  +  »/(0,  ^ 

(6.5) 

where  w  and  are  respectively  the  axial  displacemen.  and  absolute  temperature,  p  k  the  mass 
density,  A  and  p  are  the  Lame  (elasticity)  .iarameters,  e  k  the  specific  heat  and  «c  k  the  thermal 
conductivity.  The  positive  constant  ^  k  a  reference  temperature-tbe  absolute  temperature  of  a 
stress-free  reference  state  for  the  ro^.  The  nonnegative  constants  op  and  oi,  are  respectively  the 
viscoelastic  coeflScient  and  the  coefiBcient  of  thermdl  e.rpnrsiwii,  7  and  v  are  the  noise  processes,  and 
the  function  60  €  ^3(0, 1)  is  given  by 

‘"W = { S:  («■') 

Because  of  the  insulated,  or  Neumann,  boundary  conditions  in  (6.4)  00  the  temperature  dis- 
tributim,  the  open-loop  system  corresponding  to  (6.1)-(6.2)  has  a  lero  eigenvalue  for  which  the 
associated  eigenvector  consists  of  sero  dkplacement  and  velocity  and  nonsero  uniform  ten4>erature 
dktributioD.  Thk  eigenvector  k  orthogonal  (in  ^3(0, 1))  to  the  control  input  function  ^  and  k  in  the 
null  space  of  the  output  operator  corresponding  to  the  measurement  in  (6.5),  so  that  the  span  of  thk 
eigenvector  k  nnoontroUable  and  unobservable.  It  follows  that  (i)  the  only  part  of  the  temperature 
dktributioo  that  can  be  controlled  or  observed  k  the  part  that  k  orthogonal  to  unifwm  tempe^ure 
distributions;  (ii)  the  average  (over  17)  temperature  in  the  rod,  which  we  denote  by  k  neither 
stabilixable  nor  detectable;  (ui)  k  a  coitstant  function  of  t. 

Consequently,  in  the  thermoelastic  control  problem,  we  replace  the  temperature  distribution 
tf(<,i7)  with 

(6.7) 

The  state  equations,  then,  are  (6.1)-(6.5)  with  9  replaced  by  9.  The  state  space  H  has  the  structure 
in  (2.20)  with 

Ho  =  L,(0,l),  Vb  =  £f^(0,l),  (6.8) 
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Vi  =  H\Q,i)C\H2. 


/fa  =  {^  €  la(0. 1)  :  [\  dr,  =  0),  K,  =  ^‘(0, 1)  n  (6.9) 

Jo 

All  of  the  spaces  in  this  example  are  real.  We  use  the  standard  inner  product  for  Ho,  but  we  use 


4rl>dr, 


(6.10) 


for  the  inner  product  on  H^.  This  inner  product  on  Hi  is  required  to  get  the  £2  for  which  the 
semigroup  generator  in  this  example  has  the  form  in  (2.19).  For  Vq  and  Vi,  we  use  the  norms 


Iloilo  =  {f^  W?  dr,yi\  ll^ll,  =  ( ‘  dq)‘/». 


(6.11) 


We  define  the  operators  Aj  €  B(V^- ,  Vj),  j  —  0,2,  Do  €  B(Vo,  VJ)  and  £o  6  B(Vo,  Vj)  by 

{d>,Aorp)o  =  t  dt},  ^,rl>eVo,  (6.12) 

Jo  P 

{d>,Airl>)i  =  C  ^4' ri>' dr,,  ^.V>€Vj.  (6.13) 

Jo  P® 

{4>,Votf>)o=  *,^€Vo,  (6.14) 

Jo  P 

and 

(^,  £ot&)a  -  -  t  dr,,  d€Vi,  ^  €  Vo.  (6.15) 

Jo  P 

With  these  operators,  the  system  in  (6.1)  and  (6.2),  with  9  replaced  by  9,  has  the  form  in  (2J22) 
with  a  semigroup  generator  of  the  form  in  (2.19). 

From  (6.12)-^6.15),  it  fc^ows  that  we  have  all  of  the  conditions  in  Section  5,  including  the 
hypotheses  of  Theorem  5.1  (assuming  at  >  0,  else  we  would  not  have  a  thermoelastic  problem). 
Hence,  the  open-loop  thermoelastic  system  is  uniformly  exponentially  stable,  even  if  ap  =  0. 

For  the  numerical  studies  in  this  paper,  we  chose  the  parameters  m  (6.1)  and  (6.2)  for  an  alu¬ 
minum  rod  of  length  lOOin  (see  [27,  29]).  With  the  length  normalised  to  1,  the  parameters  take  the 
values  in  Table  6.1. 


p  =  9.82  X  10“’ 

A  =  2.064  X  10-* 

p  =  1.11  X  10-* 

e  =  5.40  X  10-=^ 

K  =  7.02  X  10-^ 

^  =  68 

at  =  1.29  X  10-’ 

Op  =  0 

in  =  J85 

17,  =  .486 

Thble  6.1:  Parameters  for  (6.1)-(6.6) 

The  numerical  results  in  this  paper  focus  on  the  effects  of  thermoelasUc  damping.  In  [7],  we 
presented  numerical  results  for  a  similar  example  that  included  nonsero  viscoelastic  damping  (ap  > 
0).  The  functional  gains  were  much  smoother  than  the  gains  for  the  case  with  thermoelastic  damping 
only,  and  the  approximating  functional  gains  converged  much  faster.  The  numerical  results  in 
[7]  indicate  that,  if  Voigt-Kelvin  viscoelastic  damping  is  present,  its  effect  dominates  the  effect  of 
thermoelastic  damping,  but  it  is  not  clear  whether  Voigt-Kelvin  viscoelastic  damping  is  present  at 
significant  leveb  in  conunon  metals. 
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6.2  The  Optimal  Control  Problem  and  the  Approximation  Scheme 

We  have  m  =  /  =  p  s  1  with  the  input  operators  given  by 

Bor  =  Bor  =  (-tolr,  r  €  R^,  Bj  =  =  0,  (6.16) 

P 

and  the  output  operator  given  by 

C{4>,il>,0)  =  i^.i>.0)€H^VoxHoxBi.  (6.17) 

In  the  quadratic  performance  index,  we  take  the  operator  Q  €  B{H)  to  be  given  by 

Qz  =  Q(w,  to,  )  =  (w,  ti,  0),  (6.18) 

and  we  take  A  =  1.  This  Q  penalises  the  total  nwyhanifial  energy  in  the  rod  but  does  not  penalise 
temperature  variations  from  the  constant  average  value.  The  operator  Q  €  B{H)  is  given  by  (3.10) 
with  B  =  B  given  by  (3.4)  and  (6.16).  Since  j(t)  and  v(t)  have  unit  intensities,  F  =  B  =  1. 

The  optimal  function^  control  and  estimator  gains  have  the  form  ki  =  (ki,x,ki,3,fci^)  and 
*1  =  (ii,iiii.ai^i.s)  wth  €  /fo(0,l),  4ri,a,ki,a  €  X-a(0, 1),  and  ki,3,ki4  6  C  1^(0, 1).  If 

K  and  K  are,  respectively,  the  control  and  estimator  gain  operators,  then 

Kx^  +  *  =  (^,^,8)€B,  (6.19) 

and 

Br  =  (rii.i,rii,3,rii,3)€B,  reR}.  (6.20) 


In  [7],  we  compared  two  Galerkin  approximations  for  solving  a  linear-quadratic  regulator  problem 
for  the  thermoelastic  rod  in  this  example.  One  scheme  was  a  finite  element  approxinaation  in  which 
linear  splines  were  the  basis  vectors;  in  the  other  approximation,  the  open-loop  eigenvectors  of  the 
distributed  systems  were  the  basis  vectors.  The  inodal  ^proximation  gave  fiister  ermvergence  for 
the  approximating  functional  control  gains.  In  this  paper,  we  use  the  modal  approximation  only. 

It  is  easy  to  see  that,  for  the  boundary  conditions  in  this  example,  the  eigenspaces  of  the  open-kx>p 
thermoelastic  rod  are  three-dimensional  subspaces  each  spanned  by  a  two-dimensional  subspace  o' 
the  undamped  wave  equation  and  a  one-dimensional  eigenspace  of  the  heat  equation.  The  eigenvec¬ 
tors  of  the  wave  equation  are  ^e  waves,  and  the  eigenvectors  of  the  heat  equation  are  cosine  waves. 
The  sequence  of  three-dimensional  subspaces  of  the  thermoelastic  rod  are  mutually  orthogonal  and 
complete  in  the  state  space  H.  Thus  it  is  easy  to  show  that  all  the  conditions  of  Hjrpothesis  4.4 
hold. 

The  open-loop  eigenvalues  can  be  determined  as  the  solutiou  to  the  cubic  characteristic  equations 
corresponding  to  the  three-dimensional  eigenspaces.  For  the  values  of  the  parameters  that  we  used, 
the  eigenvalues  corresponding  to  each  open-loop  subspace  consist  of  a  coii^>lex  cotyugate  pair  and 
a  real  eigenvalue,  all  with  negative  teal  parts.  It  can  be  shown  by  anabsis  of  the  sequence  of 
cubic  equations  that,  asymptotically,  the  teal  eigenvalues  approach  — oo  and  the  complex  pairs  of 
eigenvalues  approach  a  vertical  line  strictly  to  the  left  of  the  ima^ary  axis.  This,  together  arith  the 
orthogonality  and  completeness  of  the  eigenspaces,  guarantees  (4.24);  Le.,  that  the  approximating 
open-loop  semigroups  are  uniformly  exponentially  stable,  with  a  decay  rate  uniform  in  n  (the  order 
of  approximatioo,  or  number  of  nx^al  subspaces).  Hence  (4.18)  and  (4.19)  hold.  Therefore,  (4.20)-’ 
(4.30)  are  guaranteed. 

To  obtain  the  approximating  control  and  estimator  gains  shown  in  Figures  6.1  and  6.2,  we  used 
the  matrix  sign  function  method  in  [30]  to  solve  Riccati  matrix  equations  equivalent  to  the  finite 
dimensional  Riccati  operator  equations  (4.1)  and  (441),  as  discussed  in  Section  4.3.  We  used  (4.74) 
and  (4.75)  with  m  s  p  s  1  to  eoii4>ute  the  approximating  functional  control  gains  t7^(i  =  1,2,3) 
and  ^>proKimating  functional  estimator  gains  *m(**  lt2,3). 
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6.3  Numerical  Results  for  Finite-Dimensional  Compensators 

In  each  of  the  figures,  we  have  plotted  the  approximation  to  the  particular  functional  g*in  for  each  n 
between  18  and  33,  where  n  is  the  number  of  modal  aubspaces  u^.  Because  the  damping  produced 
by  thermoelastic  dissipation  is  so  small  in  this  example,  we  see  nothing  resembling  gain  convergence 
until  we  use  at  least  n  s  15.  The  convergence  results  for  approximations  to  the  infinite-dimensional 
LQG  problem  guarantee  that  all  of  the  functional  gains  do  converge,  but  the  convergence  theory  does 
not  indicate  the  rate  at  which  the  gains  converge.  Numerical  e3q>erience  has  shown  that,  generally, 
greater  damping  causes  faster  gain  convergence. 

We  ate  not  sure  that  we  are  seeing  convergence  in  Figure  6.1.  Increasing  n  past  40  does  not 
make  the  functional  control  gains  look  closer  to  any  limit,  and  between  n  =  40  and  n  =  50,  the 
numerical  solution  to  the  Riccati  equation  is  so  inaccurate  in  some  cases  that  the  corresponding 
gains  do  not  resenoble  those  in  Figure  6.1.  While  the  functional  gains  must  converge,  it  is  possible 
that  the  order  approximation  required  for  convergence  exceeds  our  capability  to  solve  the  Riccati 
equations  accurately.  Another  reason  that  we  question  whether  our  plots  of  the  functional  control 
gains  show  convergence  is  that  when  we  compute  the  control  gains  for  both  aj)  =  0  and  oi,  =  0, 
the  plots  look  identical  to  Figure  6.1.  But  with  no  damping  for  the  wave  equation  and  the  coercive 
weighting  that  we  place  on  the  solution  to  the  wave  equation  in  the  performance  index,  the  norms 
of  the  finite  dimensional  Riccati  operators  are  guaranteed  to  grow  without  bound  as  n  increases 
[31,  18].  Indeed,  when  ao  =  0  and  at  =  0,  our  numerical  solutions  to  the  Riccati  equations  break 
down  for  smaller  n  than  they  do  when  =  0  and  at  >  0.  There  is  some  difference  between  the 
finite-dimensional  gain  matrices  that  we  compute  with  and  without  thermoelastic  damping  in  the 
plant  model,  but  that  difference  is  too  small  to  be  seen  in  plots  of  the  functional  gains. 

The  question  arises,  then,  whether  the  very  light  structural  damping  produced  by  the  thennoe- 
lastic  effect  in  the  rod  is  significant  in  compensator  design.  To  address  this  question,  we  computed 
eigenvalues  for  two  closed-loop  systems.  Each  closed-loop  system  was  constructed  by  connectmg  a 
compensator  based  on  a  control  model  consisting  of  the  first  20  modal  subspaces  to  a  simulaition 
model,  or  truth  model,  consisting  of  the  first  30  modal  subspaces  of  the  rod.  Each  compensator 
thus  has  dimension  60  while  the  simulation  model  has  dimension  90.  The  30-mode  Emulation  model 
was  the  sanae  in  each  case;  it  had  the  parameters  in  Ibble  6.1,  including  at  s  1.29  x  10~^.  The 
20-mode  control  model  for  Compensator  1  also  had  the  parameters  in  Ibble  6.1.  The  control  model 
for  Con:q>en8ator  2  had  at  =  0,  and  all  of  the  other  parameters  had  the  values  in  Tbble  6.1.  This 
means  that  there  is  no  damping  for  the  mechanical  vibrations  of  the  rod  in  the  open-loop  control 
model  for  Compensator  2.  Because  the  temperature  distributi(»  is  not  penalized  in  the  performance 
index,  the  control  gains  ki^  and  and  estimator  gains  and  are  all  sero  in  Compensator 

2,  and  the  gains  ki,i  and  k"  j,  ki,i  and  k|  i,  ki,3  and  k”],  and  k”]  are  those  that  would  be 
computed  for  a  20-mode  nnodel  of  the  undamped  wave  equation  alone. 

Table  6.2  shows  typical  eigenvalues  for  the  open-loop  system  and  for  the  closed-loop  system 
produced  by  each  compensator.  Since  each  compensator  contains  a  copy  of  each  of  the  first  20 
modal  subspaces,  each  closed-lo^  system  contains  six  states,  and  six  eigenvalues,  corresponding  to 
each  of  the  first  20  modal  subspaces.  Each  closed-loop  system  also  contains  the  30  states  in  twenty- 
first  through  thirtieth  modal  subspaces.  While  the  closed-loop  performance  in  the  first  ten  or  so 
modes  is  similar  with  both  CMnpensators,  the  closed-loop  eigen^ues  corresprmding  to  several  of  the 
higher-frequency  modes  reveal  important  differences  between  the  two  omnpensators.  In  particular, 
we  note  the  second  complex  pair  of  closed-loop  eigenvalues  listed  for  mode  18.  The  magnitude  of  the 
real  part  produced  by  Compensator  1  is  more  than  20  times  the  corresponding  number  produced  by 
Compensator  2.  The  same  is  true  for  mode  19.  In  certain  high-frequency  closed-loop  states,  then, 
the  decay  rates  produced  by  Compensator  1  are  more  than  20  times  the  decay  rates  produced  by 
Compensator  2. 

The  eigenvalues  in  Ibble  6.2  for  modes  21  and  22,  the  first  modes  not  modelled  in  the  com¬ 
pensators,  are  typical  of  the  eigenvalues  for  all  ten  modes  that  are  present  in  the  simulation  model 
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but  not  in  the  control  models.  These  eigenvalues  show  that  we  have  modelled  enou^  modes  in  the 
compensators  ‘  >  eliminate  any  significant  spillover  between  modelled  and  unmodelled  modes. 

Because  the  magnitudes  of  the  real  eigenvalues,  which  correspond  to  the  heat  equation  (6.2), 
are  so  much  larger  than  the  magnitudes  of  the  complex  eigenvalues,  we  suspected  that  it  might  be 
possible  to  eliminate  the  states  corresponding  to  the  real  open-loop  eigenvalues  from  the  control 
model  and  base  a  compensator  design  on  a  control  model  crmsisting  of  a  sequence  of  secmd-order 
modes  with  eigenvalues  equal  to  the  complex  open-loop  eigenvalues  of  the  thermoelastic  rod.  This 
amounts  to  putting  artificial  viscous  damping  in  the  wave  equation. 

We  carried  out  such  a  design  with  twenty  second-order  modes  having  eigenvalues  equal  to  the  first 
twenty  pairs  of  complex  open-lo<^  thermoelastic  eigenvalues  and  mode  shapes  the  same  as  the  first 
twenty  nnodes  of  the  undamped  rod.  This  compensator  had  dimension  40.  When  we  closed  the  loop 
with  the  30-mode  simulation  model  used  for  Table  6.2  and  computed  the  closed-loop  eigenvalues,  we 
obtained  virtually  identical  results  to  those  for  Compensator  1,  except  that  this  third  closed-loop 
system  had  only  half  as  many  real  dgenvalues  because  the  corresponding  states  were  not  modelled 
in  the  compensator.  Even  for  modes  18  and  19,  all  of  the  closed-loop  eigenvalues  produced  by 
the  third  compensator  matched  to  at  least  three  digits  the  corresponding  eigenvalues  produced  by 
Compensator  1. 
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Figure  6.2a:  Approxunatmg  F^mctiooal  Ectimator  Gains  k*  i ,  n  =  18, 19, ...  33 
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Mode 

Number 

Open-loop 

Cloeed-loop 
with  Compensator  1 

Closed-loop 
with  Compensator  2 

1 

-2.30  X  10-^±i6.57x  10° 

-1.30  X  lO-* 

-3.15  X  10-*  ±»6.57x  10° 
-1.44  X  10°  ±i  6.89  x  10° 
-1.30  X  10-“ 

-1.30  X  10-“ 

-3.13  X  10-*±t6.58x  10° 
-1.45  X  10°  ±16.87  X  10° 
-1.30  X  10-“ 

-1.31  X  10-“ 

2 

-9.19  X  10-^  ±1 1.31  X  10‘ 

-5.21  X  lO--* 

-1.26  X  10-*  ±i  1.31  X  10* 
-1.22  X  10-*±i  1.31  X  10* 
-5.21  X  10-“ 

-5.21  X  10-“ 

-1.63  X  10-*±il.31  X  10* 
-8.50x  10-3  ±1 1.32  X  10* 
-5.21  X  10-“ 

-5.23  X  10-“ 

3 

-2.07  X  10-«±il.97x  10* 

-1.17  X  10-3 

-2.55X  10-*  ±11.97  X  10* 
-3.45x  10-*±*1.97x  10* 
-1.17  X  10-3 
-1.17  X  10-3 

-2.17  X  10-*±il.97x  10* 
-3.84  X  10-*±il.96x  10* 
-1.18  X  10-3 
-1.17  X  10-3 

10 

-2.30  X  10"*  ±16.57  X  10* 

-1.30  X  10-3 

-1.82  X  10-*  ±  16.57  X  10* 
-8.03  X  10-3  ±»  6.57  X  10* 
-1.30  X  10-3 
-1.30  X  10-3 

-2.16  X  10->  ±i6.56  x  10* 
-4.66  X  10-3  ±i 6.58  X  10* 
-1.30  X  10-3 
-1.31  X  10-3 

14 

-4.50  X  10-®±i9.20  x  10* 

-2.55  X  10-3 

-3.44  X  10-3  ±i 9.20  X  10* 
-3.41  X  10-3  ±t  9.20  X  10* 
-2.55  X  10-3 
-2.55  X  10-3 

-3.76  X  10-3  ±i  9.19  X  10* 
-2.66  X  10-“  ±i 9.20  x  10* 
-2.56  X  10-3 
-2.55  X  10-3 

18 

-7.45  X  10-*±t  1.18  X  103 

-4.22  X  10-3 

-1.53  X  10-3  ±1 1.18  X  103 
-2.16  X  10-3  ±i  1.18  X  103 
-4.22  X  10-3 
-4.22  X  10-3 

-1.74  X  10-3  ±i  1.18  X  103 
-9.75x  10-*±il.l8x  10* 
-4.22  X  10-3 
-4.24  X  10-3 

19 

-8.30  X  10-*±»1.25x  103 

-4.70  X  10-3 

-1.03  X  10-3  ±i  1.25  X  103 
-1.98  X  10-3  ±»  1.25  X  103 
-4.70  X  10-3 
-4.70  X  10-3 

-1.22  X  10-3  ±i  1.25  X  io» 
-9.05  X  10-®±il.25  X  103 
-4.70  X  10-3 
-4.72  X  10-3 

20 

-9.19  X  10-*±il.31  X  103 

-5.21  X  10-3 

-2.51  X  10-3  ±«  1.31  X  10’ 
-6.67X  10-“±il.31  X  103 
-5.21  X  10-3 
-5.21  X  10-3 

-3.18  X  10-3  ±1 1.31  X  10* 
-9.33  X  10-*±il.31  X  103 
-5.21  X  10-3 
-5.23  X  10-3 

21 

-1.01  X  10-“*  ±1 1.38  X  10* 
-5.75  X  10-3 

-1.04X  10-“±il.38x  103 
-5.75  X  10-3 

-1.04  X  10-“±iU8  X  103 
-5.75  X  10-3 

22 

-1.11  X  10-'‘±il.45x  103 
-6.31  X  10-3 

-1.34  X  10-“±il.45x  103 
-6.31  X  10-3 

-1.34  X  10-“  ±i  1.45  X  103 
-6.31  X  10-3 

Simulation  Model: 
Compensator  1: 
Compensator  2: 

30  Modal  Subspaces,  =  1.29  x  10-3 

20  Modal  Subspaces,  at  =  1.29  x  10-3 

20  Modal  Subspaces,  at  =  0 

Table  6.2:  Typical  Open-loop  and  Cloeed-loop  Eigenvalues 
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7  Conclusions 


The  abstract  formulation  of  distributed  models  and  the  approximation  theory  developed  in  this 
paper  apply  to  a  wide  variety  of  thermoelastic  control  systems.  The  uniform  exponential  stability 
result  in  Section  5  applies  to  a  large  class  of  thermoelastic  problems,  but  not  to  certain  systems  that 
are  known  to  be  uniformly  exponentially  stable  [14]. 

The  numerical  study  in  Section  6  focussed  on  the  effect  of  thermoelastic  damping  in  optimal 
control  of  a  flexible  structure.  The  eigenvalue  results  demonstrate  that,  even  though  thermoelastic 
damping  is  small  in  common  metals,  a  compensator  baaed  on  a  thermoelastic  model  of  a  flexible 
structure  can  produce  signiflcactly  better  response  in  high-frequency  modes  than  a  ccxnpensator 
based  on  an  undamped  model  can  produce. 

The  theory  in  Sections  2-4  also  applies  to  thermoelastic  control  problems  in  which  a  thermal 
disturbance  excites  mechanical  vibrations.  This  class  of  problems,  which  includes  vibrations  in 
flexible  space  structures  caused  by  solar  heating,  might  provide  the  most  important  applications  for 
the  theory  developed  here. 
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This  paper  derives  sufficient  conditions  for  nniforai  exponentiaJ  stability  of  solutions  to 
abstract  linear  evolution  equations  that  are  second-order  in  time.  The  main  problems  motivating 
the  paper  involve  sets  of  coupled  partial  differential  equations  with  nonsymmetric  damping.  Such 
systems  arise  in  acoustics  problems  when  the  interaction  between  compressible  fluids  and  elastic 
boundaries  produces  a  skew-symmetric  damping  term. 
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1  Introduction 


By  abstract  wave  equations,  we  mean  evolution  equations  that  are  second-order  in  time.  This 
class  includes  both  classical  wave  and  biharmonic  equations.  The  evolution  equations  that  are  the 
primary  motivation  for  the  paper  are  sets  of  coupled  partial  differential  equations  that  represent 
acoijstic  waves  in  a  compressible  fluid  interacting  with  an  elastic  boundary. 

This  paper  has  two  main  objectives;  first,  to  derive  sufScient  conditions  for  the  solutions  to 
abstract  linear  wave  equations  with  coercive  damping  to  be  uniformly  exponentially  stable;  second, 
to  demonstrate  that  the  class  of  coupled  partial  differential  equations  of  greatest  interest  in  this 
paper  do  indeed  have  the  abstract  second-order  form  treated  here.  Uniform  exponential  stability  is 
important  in  active  control  of  distributed  systems  [1,2,3],  and  having  the  abstract  second-order  form 
allows  an  extensive  theory  of  approximation  derived  especially  for  control  and  estimation  problems 
to  be  applied  (see  [1,2,3,4],  for  example). 

The  most  interesting  property  of  the  equations  motivating  this  paper  is  that  the  damping  op¬ 
erator  is  not  symmetric.  A  nonzero  skew-symmetric  damping  operator  makes  the  issue  of  uniform 
exponential  stability  more  complicated  than  in  the  case  of  s>‘mmetric  damping.  In  particular,  when 
the  damping  operator  is  symmetric,  it  must  only  be  coercive  with  respect  to  the  kinetic-energy  norm 
to  guarantee  uniform  exponential  stability.  However,  an  example  in  Section  3  of  this  paper  shows 
that,  in  general,  more  is  needed  when  the  damping  is  not  symmetric.  (While  the  operator  in  the 
term  with  the  first-order  time  derivative  in  second-order  evolution  equations  is  referred  to  commonly 
as  damping,  the  skew-symmetric  part  of  this  operator  represents  conservative  energy  transfer  among 
different  states,  as  opposed  to  dissipation.) 
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2  Abstract  Linear  Wave  Equations 

Throughout  this  paper,  H  or  Hj  (j  =  0, 1,2)  will  be  a  Hilbert  space  with  inner  pro'luct  (■,  •)  or  {-,  ■}j 
and  corresponding  induced  norm  |  ■  1  or  |  ■  [,  .  Also,  V  or  V)  will  be  a  Hilbert  space  with  norm  ||  ■  |{ 
or  II  •  llj  and  inner  product  {■,  )v  or  (  ,  The  continuous  dual  of  V  will  be  denoted  by  V',  and 

V^H^V  (2.1) 

will  mean  that  V  is  embedded  densely  and  continuously  in  H,  which  implies  that  H  is  embedded 
densely  and  continuously  in  V.  In  this  case,  (•,  )  will  denote  both  the  H-inner  product  and  the 
duality  pairing  on  V  x  V'. 

For  an  operator  C  G  B{Vi,VJ),  we  define  £*  G  B(Vj,V/)  (the  Banach-space  adjoint  of  £)  by 

(ii.CV),  =  (u;,£v)^,  v€Vi,w£Vj.  (2.2) 

We  say  that  C  G  B(V,  V’)  is  symmetric  if  £  =  £*.  We  say  that  £  G  B(V,  V)  is  nonnegative  if 

Re(v,£v)>0,  vEV,  (2.3) 

and  that  £  is  dissipative  if  — £  is  nonnegative.  We  say  that  £  G  B{V,  V)  is  H-coercive  if  there  exists 
a  positive  real  number  p  such  that 

.^e(i',  £v)  >  ;i|t;p,  V  £  V,  (2.4) 

and  that  £  G  B{V,  V)  is  r-coen;t-e  if  there  exists  a  positive  real  number  n  such  that 

Re(v,£t;)>^M|^  f  G  K  (2.5) 

We  assume 

Vo  -  Ho  Vo'.  (2.6) 

and  we  study  evolution  equations  of  the  form 

U'(i)  +  r>od)(0  +  •4otn(t)  =  0,  1  >  0,  (2.7) 

where  Vo,Ao  €  B{\’o,  V^]  with  Po  nonnegative  and  ^o  equal  to  the  Riesz  map  on  Vq.  Thus, 

{v,w)v„  =  {v,Aow)o,  (2.8) 

and  -4o  is  symmetric  and  t  'o-coercive 
We  define 

V  =  Vo  X  Vo  --  H  =  Vo  X  Ho  --  V'  =  Vo  X  VJ.  (2.9) 

(We  identify  Vq  with  VJ  in  the  first  component  of  H  only.)  We  define 

If  follows  from 

]  €B(V'/.V,)  (2.11) 

that  v4  is  an  isomorphism  from  V  to  V'.  Now  we  define  A  :  Dom(A)  C  H  — *  H  by 

Dom(.4)=  {(v,/.)G  V  :>((v,/i)gH},  >4  =  ^lDon,(yt)-  (2.12) 

According  to  Lemma  2.1  in  [3],  A  is  a  maximal  dissipative  operator  on  H.  Hence  A  generates  a 
strongly  continuous  contraction  semigroup  on  H.  (See  [4,2,5]  for  Bimilar  approaches.) 
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3  Uniform  Exponential  Stability 

Theorem  3.1  If  T>o  is  Ho-coercive  and  there  exists  a  positive  real  number  7  such  that 

IR^^.PoMol  <  tIIHIo  •  lR^<A.:DoA)oP'^  v.heVo,  (3.1) 

then  the  semigroup  generated  on  H  by  the  operator  A  in  (2.12)  is  uniformly  exponentially  stable. 
Proof  We  define  the  following  self-adjoint  bounded  linear  operator  on  II: 

<3  =  [  "/  ]  (3-2) 

where  a  is  a  positive  real  number.  For  e  sufficiently  large,  Q  is  /f-coercive.  Also,  QV  C  V.  For 
X  =  (t),  h)  €  Dom(y4)  C  V, 

Re  {Qx,  Ax)  =  Re  {Qx,Ax)  =  -||t;llo  +  l^lo  “  <7Re(h,Pof»)o  -  Ra{v,'Doh)o.  (3.3) 

Since  Vq  is  Ho-coercive,  it  follows  from  (3.1)  and  (3  3)  that,  for  cr  suificicntly  large,  there  exists  a 
positive  real  number  p  such  that 

Re{Qx,  Ax)  <  — /i|x|^,  x  €  Dom(i4).  O  (3.4) 


Remark  3.2  A  sufficient  condition  for  the  existence  of  a  positive  real  number  -)  such  that  (3.1) 
holds  IS  that  Vo  be  Vo-coercne. 

Remark  3.3  If  Vo  is  Ho-coercive,  the  generalized  Schuarz  inequality  and  the  fact  that  Vo  € 
B(Vo,Vo)  imply  that  |Re(i  , (Po  +  Po)^)ol  *s  bounded  by  the  right  side  o/ (3.1)  for  some  positive 
7  independent  of  t  and  h.  Therefore,  when  Vq  w  Ho-coercive,  (3.1)  holds  for  some  7  independent 
of  V  and  h  if  and  only  if 

|Re(v,(Po-Po)h)o|<7iMlo■|Re{h,Pch)oI''^  v,h€Vo,  (3.5) 

for  some  7  independent  of  v  and  h.  When  Vo  is  Ho-coerctve,  one  sufficient  condition  for  (3.5)  is 
that  7l{Vo  —  Vq)  C  Ho  (this  includes  the  case  Vo  —  Vq). 

The  following  example  shows  that  Vo  being  ffo-coercive  without  (3.1)  or  some  other  condition 
is  not  sufficient  for  the  semigroup  generated  by  A  to  be  uniformly  exponentially  stable. 

Let  a  be  a  negative  real  number,  and  let  and  Wn  be  sequences  of  positive  real  numbers  such 
that 

u„  —•  00  and  0nfio„  —•00  as  n  — •  00  (3.6) 

and 

8up^„/w’<oo.  (3.7) 

n 

Let  Ho  be  tj,  let  Tq  be  the  infinite-dimensional  diagonal  matrix  with  the  sequence  Wn  on  the  diagonal, 
let  Vq  =  T^^Ho  ^ith  Ao  =  Tq,  and  let  Vq  be  the  infinite-dimensional  diagonal  matrix  with  the 
sequence  a  -f  i0n  on  the  diagonal. 

The  eigenvalues  of  the  corresponding  semigroup  generator  A  are  the  solutions  to  a  sequence  of 
quadratic  equations.  Elementary  analysis  shows  that  these  eigenvalues  approach  the  imaginary  axis 
as  n  — •  00.  Hence  the  semigroup  generated  on  /f  by  A  is  not  unif(»^mly  exponentially  stable. 
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4  Coupled  Abstract  Wave  Equations 

In  this  section,  we  assunie  that  the  Hilbert  spaces  Ho  and  Vb  and  the  operators  Ao  and  Vo  have  the 
forms 


Ho  =  Hix  Hi, 

Vo  =  1 

X  Vi, 

(4.1) 

Ao  = 

f  ^11 

1  0 

0 

An 

t 

(4.2) 

Vn 

-vii 

Vn 

Vn 

]■ 

(4.3) 

We  assume  that  An  is  the  Riesz  map  for  t  =  1,2,  and  that  1>,^  e  with  Vn  nonnegative. 

We  are  particularly  interested  in  problems  where  the  wave  equation  on  Vi  x  Hi  represents  an 
elastic  boundary  condition  for  the  wave  equation  on  x  H^.  In  this  case,  we  have 

HiVii)  C  Hi.  (4.4) 

We  should  note  that  (4  4)  does  not  imply 

7^(DI,)  C  Hi.  (4.5) 

We  consider  the  domain  of  the  semigroup  generator  A  for  the  coupled  system  when  (4.4)  holds. 
It  follows  from  (2.9)-(2.12)  that  Dom(.4)  is  the  set  of  all  elements 

ivuVi.,hi,hi)€V  =  Vi  X  Vjx  Vi  X  Vj  (4.6) 

such  that 

AiiVi -i-Viihi  ^  Hi,  (4.7) 

AiiVi +  V2ihi —  Vl^hi  ^  Hi.  (4.8) 

Frc.n  (4.7),  we  see  that  the  conditions  on  vi  and  hi  are  independent  of  the  coupling  operator  Vn 
and  are  the  sane  as  when  the  two  wave  equations  are  uncoupled.  On  the  other  hand,  if  Tl{Dn)  “ 
not  contained  in  Hi,  then  the  conditions  on  Vi  and  /tj  do  depend  on  Vn.  In  most  applications.  PJj 
affects  the  ‘natural  boundary  conditions’  for  Vi  and  hi. 

To  be  more  concrete,  we  take  Tj  €  B{Vi,Hi)  such  that 

{\.w)v:,=-{v.AiiM:)i  =  {TiV,Tiw)i,  v,u;6V'2,  (4.9) 

and  we  take  7^  to  the  Hilbert  space  adjoint  of  Ti  with  respect  to  the  /fj-inner  product.  It  is  a 
straightforward  exercise  lo  show  that  (4.6)  is  equivalent  to 

TiVi  +  TiA^iVahi  ~  TiAj^'^nhi  €  Dom(T2).  (4.10) 

In  many  examples,  it  is  easy  to  determine  TiA^iVn  and  TiA^Vn-  If  'R.{Vn)  C  Hi,  the  term 
TiA^iViihi  drops  out  of  (4.10),  but  this  ts  not  the  case  in  which  we  will  be  interested  from  here 
on. 

Before  considering  an  example,  we  will  give  a  sufficient  condition  for  solutions  to  the  system  of 
coupled  abstract  wave  equations  to  be  uniformly  exponentially  stable  ia  H  =  Vo  x  Ho. 

Theorem  4.1  Let  Ho,  Vo,  Ao  and  Vo  have  the  forms  in  (4.1)-(4.3)  mih  Vn  symmetric,  i  =  1,2. 
Lei  (4.4)  hold.  IfVii  is  Hi-coercive  and  Vn  is  V]~coercive,  then  the  semigroup  generated  on  H  by 
the  operator  A  in  (2.9)-(2.12)  is  uniformly  exponentially  stable. 

Proof  Since  Vn  €  B(Vj.  Vj),  (4.4)  implies  Vn  €  B(Vj,  Hi).  For  wj, hi  €  Vi  and  Vi,hi  6  Vj. 

+  {^,'Pnhi)i\  =  +  (hx,I>ijV3)il 

<  ids'll!  •  IIAjIIj  +  \hi\i  -  ||wj||j)-  (4.11) 

Since  Vn  is  svinmetric  and  //j-coercive,  Vn  is  symmetric  and  Vj-coercive  and  the  Vj-norm  is 
stronger  than  the  //i-norm,  it  follows  from  (4.11)  that  (3  5)  bt^ds  for  some  postive  7  independent 
of  u  =  (vj ,  Vi)  €  Vq  and  />  =  (bj ,  hi)  €  O 
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5  An  Example  from  Acoustics 

In  this  one-dimeosioDal  application,  the  wave  equation  on  Vj  x  Hi  governs  the  velocity  potential 
in  a  compressible  fluid,  and  the  wave  equation  on  Vi  x  Hi  reduces  to  the  equation  of  motion  of  a 
mass-spring-damper  system  on  one  boundary  of  the  fluid.  We  assume  zero  velocity  potential  on  the 
other  boundary  of  the  fluid. 

We  have 


Hi  =  Vi=R\ 

Hi  =  1,(0, 1),  Vi  =  {<f>e  H\Q,  1) :  ^(1)  =  0), 

-4ii  =  1 

=  f  dr) 

Jo 

=  <1. 

Vi2  =  C2.4„, 

{0,'DiiiP)i  =  -^,^.(0),  0eVi,xt>€  V,. 

We  have  taken  all  physical  constants  to  be  1,  except  the  nonnegative  real  numbers  fj  and  f,. 
The  operator  T,  in  (4.9)  and  (4.10)  is  given  by 

Dom(Tj)  =  Vi,  Ti4>  =  4>' , 

and  its  adjoint  is  given  by 

Dom(7^)  =  (d  €  /f‘(0,l)  :  d>(0)  =  0}. 

In  this  example,  is  a  function  in  £,(0,1).  From  (5.7),  it  follows  that 


(5.1) 

(5.2) 

(5.3) 

(5.4) 

(5.5) 

(5.6) 

(5.7) 


=  0,  0<i7<l,  0£RK 


Also, 


(5.8) 

(5.9) 

(5.10) 

AiiVii  =  dl.  (5.11) 

It  follows  from  (4.6),  (4.7)  and  (4.10),  then,  that  Dom(y4)  is  the  set  of  all  quadruples  {q,o,0,i^) 
satisfying 

a,0eR\  (5.12) 

o.vl.e/f‘(0,l),  <>(!)  =  v(l)  =  0,  (5.13) 

{<^  +  (irP-0y  €H\0.\),  <>'(0)  +  f  jV-HO)  -  ^  =  0.  (5.14) 

From  (2.10)-(2.12)  and  the  deflnitions  of  the  various  operators  in  this  example,  it  follows  that,  for 
{q,<I>,0,\(>)  €  Dom(y4), 


a  \  /  0 

^  1  _  ^ 

0  I  -Q  -  <1^ -f- tl>(0) 

\  ^  /  \  (^'  +  <J^')' 


(5.15) 


We  define  iy(<)  =  (a((),^(<))  €  Vo  =  Vj  x  Vj  to  be  the  solution  to 


t  >0. 


(5.16) 
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Equivalently,  a(t)  and  ^(t)  =  </>{(,  v)  satisfy 


o(t)  +  fid(0  +  Q(t)  =  ^{t,  0) 

t  >0, 

(5.17) 

d^4>,  .  d  .86  8^6 

0<i7<l,  f>0. 

(5.18) 

t>o, 

(5.19) 

^(<,1)  =  0,  t>0. 

(5.20) 

If  <1  and  fj  are  positive,  then  Theorem  4.1  says  that  the  semigroup  generated  on  ^  by  A  and 
consequently  all  solutions  to  (5.17)-(5.20)  are  uniformly  exponentially  stable.  It  is  important  to 
note  that,  to  apply  Theorem  4.1,  we  need  the  V27  in  (5.6),  or  something  similar,  as  oppKxed  to,  say, 
V22  =  (2/.  The  V22  in  (5.6)  is  a  realistic  dissipation  term  for  waves  in  a  compressible  fluid  because 
it  represents  viscosity  [6,  page  315],  (7). 

6  Conclusions 

It  is  straightforward  to  generalize  the  example  here  to  wave  equations  in  higher  dimensions.  For 
example,  the  wave  equation  on  Vi  x  Ni  could  represent  an  elastic  membrane  interacting  with  the 
waves  in  a  three-dimensional  volume  of  fluid,  represented  by  the  wave  equation  on  x  ffj.  In  such 
cases,  the  product  on  the  right  side  of  (5  7)  becomes  an  1.2-iQDer  product  on  the  elastic  boundary. 

In  the  higher-dimension  problems,  the  damping  operator  Vn  for  the  elastic  boundary  need  be 
only  /fi-coercive  for  uniform  exponential  stability,  if  V22  represents  viscosity,  as  in  the  example  here. 
This  is  important  because  many  common  dampmg  models  for  flexible  structures  are  i^i-coercive  but 
not  Vi-coercive. 

Whether  the  hypotheses  in  Theorem  3.1  or  the  requirement  in  Theorem  4.1  that  P22  be  V2- 
coercive  can  be  weakened  is  an  open  and  interesting  question.  The  example  at  the  end  of  Section 
3  suggests  that  significantly  weaker  sufficient  conditions  than  those  given  in  this  paper  for  uniform 
exponential  stability  might  not  exist. 
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ABSTRACT 

Adaptive  lattice  algorithms  are  derived  for  solution  of  unwindowed  least-squares  estimation  problems 
for  AR  and  FIR  models.  The  basic  approach  is  to  embed  the  unwindowed  problem  in  a  larger  prewin¬ 
dowed  problem  and  then  eliminate  superfluous  terms  in  the  lattice.  Initializations  are  given  to  allow  the 
lattice  to  use  no  initial  parameter  estimates  or  to  include  initial  parameter  estimates  with  desL'ed  'veightings 
in  the  quadratic  criterion  for  parameter  estimation.  A  numerical  example  is  given. 
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1.  Introduction 


This  paper  treats  two  initialization  problems  for  unnormalized  lattice  implementation  of  recursive 
least-squares  estimation  of  AR  and  FIR  models.  The  first  problem,  usually  called  unwindowed  estimation, 
is  to  obtain  the  exact  least-squares  fit  to  data  when  the  data  preceding  the  initial  point  used  cannot  be  as¬ 
sumed  to  be  zero.  The  second  problem  is  to  penalize  deviation  from  initial  parameter  estimates.  Although 
solutions  to  both  problems  are  straight-forward  in  the  classical  recursive  least-squares  algorithm,  solutions 
via  lattices  require  not  only  special  initializations  but  more  complicated  lattices  than  the  conunon  prewin¬ 
dowed  lattice  [1,2,3, 4]. 

The  terms  prewindowed  and  unwindowed  refer  to  how  data  is  handled  in  estimation,  rather  than  to  the 
AR  or  FIR  model.  In  prewindowed  estimation,  all  data  before  some  initial  time  is  assumed  to  be  zero. 
Although  this  assumption  usually  is  not  correct,  it  affects  only  the  first  few  terms  in  the  cumulative  least- 
squares  criterion,  and  therefore  the  affect  on  estimated  parameters  and  predicted  data  fades  as  the  number 
of  processed  data  points  increases.  Whether  the  error  caused  by  the  prewindowing  assumption,  when  it  is 
incorrect,  is  tolerable  depends  on  the  application.  This  error  is  eliminated  by  unwindowed  estimation, 
where  the  criterion  to  be  minimized  tries  to  fit  only  true  data  with  the  AR  or  FIR  model.  In  unwindowed 
estimation,  no  assumption  is  made  about  data  preceding  that  used. 

Including  initial  parameter  estimates  makes  practical  sense  only  in  unwindowed  estimation.  Initial  pa¬ 
rameter  estimates  can  be  included  in  prewindowed  estimation,  but  their  effect  and  that  of  the  error  induced 
o :  the  prewindowing  assumption  fade  at  the  same  rate.  Also,  to  include  initial  parameter  estimates  in  a 
lattice  filter  for  an  AR  model,  the  model  must  be  embedded  in  an  FIR  model.  This  is  necessary  to  preserve 
the  shift  structure  of  the  recursion  vectors  for  the  AR  lattice. 

Both  problems  addressed  in  this  paper  have  been  addressed  with  previous  lattice  filters.  Normalized 
lattices  for  unwindowed  estimation  were  given  in  [S,6].  An  urmormalized  lattice  for  the  scalar-channel 
(single  experiment)  unwindowed,  or  covariance  problem  was  given  in  [2].  An  urmormalized  lattice  and 
initializations  for  unwindowed  RLS  estimation  with  and  without  initial  parameter  estimates  was  proposed 
in  [7],  although,  as  discussed  in  Section  5  of  this  paper,  we  have  concluded  that  the  lattice  and  initializations 
in  [7]  do  not  solve  unwindowed  problems.  While  the  lattice  here  and  the  initialization  of  parameter  esti¬ 
mates  are  different  from  those  in  [7],  the  basic  idea  of  using  the  initial  parameter  estimates  as  initial  data 
for  the  output  charmel  of  an  FIR  model  is  common  to  Section  4  of  this  paper  and  [7]. 
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The  basic  idea  behind  the  unwindowed  AR  lattice  in  Section  3  of  this  paper  is  to  use  an  artificial 
measurement  channel  in  the  prewindowed  lattice  in  [3].  This  idea  is  straightforward,  and  it  has  been  ex¬ 
plored  previously  (see  [8]).  However,  the  additional  channel  makes  the  lattice  more  complex.  The  most 
important  objective  in  developing  the  AR  lattice  is  to  exploit  the  special  structure  associated  with  the  arti¬ 
ficial  channel  to  reduce  the  expanded  lattice  to  an  efficient  solution  of  the  unwindowed  least-squares  esti¬ 
mation  problem.  In  particular,  expanding  the  lattice  with  the  artificial  channel  increases  the  dimension  of 
the  coefficient  matrices  that  must  be  inverted  in  the  prewindowed  lattice.  By  using  the  structure  induced 
by  the  artificial  channel  and  generalizing  a  well-known  formula  for  low-rank  updates  of  matrix  inverses, 
we  obtain  an  unwindowed  lattice  in  which  the  matrices  to  be  inverted  have  the  same  dimension  as  the 
corresponding  matrices  in  the  prewindowed  lattice. 

The  authors  of  [8]  concluded  that  the  normalized  lattice  they  derived  using  an  artificial  channel  was  not 
competitive  in  terms  of  computational  requirements  with  the  covariance  lattice  in  [5],  whereas  we  believe 
that  the  lattice  developed  here  will  be  comp>etitive  with  any  unnormalized  lattice  for  solving  the  unwindowed 
problem.  These  different  conclusions  about  the  two  lattices  based  on  the  artificial -channel  idea  appear  to 
stem  in  large  part  from  the  ability  to  reduce  the  dimension  of  the  matrices  that  must  be  inverted  in  our 
lattice  and  the  inability  to  reduce  the  dimension  of  the  matrices  of  which  square  roots  must  be  computed 
in  the  lattice  in  [8]. 

The  remainder  of  the  paper  :s  organized  as  follows.  Section  2  defines  unwindowed  and  prewindowed 
estimation  problems  for  an  .AR  model,  embeds  the  unwindowed  problem  in  the  prew  indowed  problem  and 
reduces  the  prewindowed  lattice  to  produce  an  efficient  solution  of  the  unwindowed  problem.  The  basic 
algorithm  for  the  entire  paper  is  Algorithm  2.1,  the  residual-error  lattice  for  the  AR  problem.  .Algorithm 
2.2  is  used  to  generate  the  estimated  AR  coefficients.  Section  3  defines  unwindowed  and  prewindowed  es¬ 
timation  problems  for  an  FIR  model  and  develops  the  additional  update  equations  that  must  be  appended 
to  the  algorithms  in  Seaion  2.  Section  4  shows  how  to  include  initial  parameter  estimates  in  the  FIR 
problem.  For  this,  the  lattice  algorithms  in  Sections  2  and  3  arc  used  but  the  time-initializations  are  replaced 
with  those  in  Section  4. 

Section  5  presents  an  e.xample  to  demonstrate  the  faster  convergence  of  the  unwindowed  lattice  versus 
the  prewindowed  lattice  and  the  advantage  of  including  initial  parameter  estimates  in  the  presence  of 
measurement  noise.  Section  6  presents  our  conclusions,  and  the  Appendix  contains  the  matrix  inversion 
lemma  that  we  use  in  Section  2. 
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2.  The  AR  Problem 

VVe  consider  an  mx  p  measurement  matrix 

AD  =  -  Aid  .  (2.1) 

where  t  is  any  integer  and  the  column  yii)  is  a  real  /n- vector,  referred  to  as  the  /“*  measurement  channel. 
The  k'^  element  of  >'‘(r)  is  called  the  k*  measurement  in  the  i"*  channel.  By  an  -order  AR  model,  we 
mean 

n 

AO  =  ^At-j)Aj  +  ^n(0,  (2.2) 

where  the  A^/s  an  px  p  matrices  referred  to  as  AR  coefficients.  We  denote  the  /“  column  of  A^j  by  A'^^. 
As  in  [3],  we  define  the  Hilbert  space 

[2=  xj  .  each  Xj^  R"'  and  ]  =  <r.  r>  <  oo)  ,  (2.3) 

where 

=  (2.4) 

/=! 

and  A  is  a  positive  real  number  called  the  forgetting  factor.  Throughout  the  paper,  <. ,  .)  means  the  inner 
product  in  (2.4).  For  ^  =  0,  1,  2,  ....  is  the  subspace  of /jiR”,/.)  such  that  each  element  of  has 
the  form  of  z  in  (2.3)  with  the  m-vector  =  0  for  j  >  k,  and 

m 

Qi^  =  orthogonal  projection  of  ,  /.)  onto  .  (2.5) 

We  denote  the  infinite  history  of  j<r)  by 

m  =  iAoAo  ■■■  =  {y\0y\t- \)y\t-2) ...  f  (2.6) 

so  that  ^\t),  the  /**  column  of  i>(0  ,  is  the  infinite  history  of  >^f),  and  utf)e  ifRf,  A).  For  any  integer  t, 
we  define  Hj,t)  =  (0),  which  is  the  zero  subspace  of /fR",  A),  and 

HAO  =  [At-  1)  --  At-  1)  At -2)  ...  At -2)  ^2.7) 

...  n)  ...  A{t  —  't)} .  =  1.  2 . 

For  A  =  0,  1,  2,  ... ,  we  define  H^t)  =  QJliO  ■  Hence  Hj^t)  =  =  H^t)  =  (0).  We  use  the 

projeaion  operators 
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P„(t)  =  orihogonal  projection  o"'  ^ -A R'^ .  /.)  onto 


(2.8) 


Pnil^)  =  orthogonal  projection  of  onto  U„jJit)  ■  (2.9) 

(We  note  PjtjOit)  =  [Pi/)i;»'(/)  -  PJ  i)  0^01  3.nd  P^  ft)  il/U)  =  [/’„.*(/)  i/'V)  -  P.JO^'U)'])  Computing 
P„_fi)u(t)  is  an  unwindoued  least-squares  problem,  and  computing  Pft)^{t)  when  \j{t)  has  a  finite  number 
of  nonzero  entries  is  a  prewindowed  least-squares  problem. 

Throughout  this  paper.  }it)  will  denote  the  mx  p  measurement  matrix  for  unwindowed  problems  and 
^(/)  will  denote  the  m  x  ^  measurement  matrix  for  prewindowed  problems.  The  histories  of  j</)  and  y(t) 
are  ipit)  and  \i/(t).  respectively.  The  projections  defined  by  (2.8)  and  (2.9)  corresponding  to><t)  and  4/{t)  are 
Pft)  and  P„  ft):  the  projection  defined  by  (2.8)  corresponding  to  j(f)  and  \}/{t)  is  Pjt). 

Problem  2.1.  (L'nwindowed  .AR  Problem)  For  t  =  1,  2,  ....  and  /?  =  0,  1.  2,  ,  t—  1,  compute 

P„..ft)\lt(t)  =  P,..-,(f)  0...  0(t);  equi'.alently.  (for  /?=  1,  2 . /—  1)  find  matrices  .-i^ft)  such  that,  for 

each((t=  1.  ...  ,  p),  the columns  of  the  matrices  minimize 

r  r 

^'•”(0  =  Y.  1^  .  (2.10) 

:=n+l 

where  |  •  |  is  the  Euclidean  norm  of  an  tn  -vector. 

Problem  2.2.  (Prewindowed  AR  Problem)  .Assume  }(t)  =  0  for  t<0  .  For  f  =  0,  1.  2,  ...,  and 

n  =  0,  1,2 . t.  compute  Pft)  it/iti:  equivalently,  (  for  n=  1,  2,  ...  ,  f)  find  matrices  y4^/t)  such  that, 

for  each  /  (/  =  1.  ...  .pi.  the  /'*  columns  of  these  matrices  minimize 

/  r 

yfV)  =  Y .^4)  -  Y^^' ■  (2.11) 

-.=0  /=; 

We  note  that,  in  Problem  2.1,  nonzero  data  for  /  ^  0  is  allowed  but  not  used,  whereas  in  Problem  2.2,  all 
data  is  assumed  zero  for  t  <  0.  The  initial  times  in  Problem  2. 1  and  2.2  are,  of  course,  arbitrar>’.  We  denote 
them  differently  so  that  we  can  embed  Problem  2.1  in  Problem  2.2  most  conveniently. 

The  lattice  in  [<SJGI.]  solves  Problem  2.2.  We  will  need  the  following  definitions  from  [&JGI.].  For 
any  integer  t  and  any  nonnegative  integer  n  , 


W  =  Ibnit)  -  ^(0]  =  [/-  P„(t+  mh-n),  (2.13) 

e^(;)  =  top  m  rows  of  f„{t) ,  forward  residual  error,  (2.14) 

r^(t)  =  top  m  rows  of  b„{t) ,  backward  residual  error,  (2.15) 

^„+i(0  =  P  P  matrix  whose  (/,  f)  element  is  (f^(t),  biit—  1))  ,  (2.16) 

^'(0  =  p  X  p  matrix  whose  (/,  y)  element  is  (fj(0,  fn(0}  .  (2.17) 

R'(t)  =  p  y.  p  matrix  whose  (/,  y)  element  is  <6^/),  6^0>  •  (218) 


The  matrices  e„(f),  r„(t),  KJ^t) ,  P'(/)  and  R’„{t)  are  used  in  the  lattice  in  [3].  The  maXnces  fj,t)  and  bj^t),  whose 
columns  are  infinitely  long,  are  used  in  deriving  the  lattice  in  [3]  but  not  in  the  fmal  algorithm.  These  in¬ 


finite  dimensional  matrices  are  used  to  derive  the  lattices  in  this  paper  also. 

To  embed  Problem  2. 1  in  Problem  2.2,  we  define 

^(0  =0,  t  <  0.  (2.19) 

MO)  =  [  0  I  /  ]  ,  (2.20) 

mxp  m^m 

MO  =  [  ><0  10  ],  t  =  1,  2 .  (2.21) 

mxp  mxm 

Thus  p  —  p  ->r  m.  It  is  easy  to  verify 

Q,-nPni^)m  0],  /=1.2.  ...  .  n^  \,2 .  t^n.  (2.22) 


This  shows  that  the  first  (f  —  n)m  rows  of  the  first  p  columns  of  PS,t)^(t)  ait  equal  to  the  first  (/  -  n)m  rows 
of  P„,r.J,t)^(t).  Since  all  rows  of  past  the  first  {t—  n)m  are  zero,  (2.22)  shows  that  the  solution 

to  Problem  2.1  for  the  sequence  y<0  is  embedded  in  the  solution  to  Problem  2.2  for  the  sequence  MO- 
Therefore,  we  can  solve  Problem  Zl  with  the  lattice  in  [3]. 

To  make  this  solution  efficient,  though,  we  must  exploit  the  special  structure  of  the  prewindowed 
problem  for>'(0-  From  (2.19)-(2.21),  it  follows  that,  for  r  =  1,  2,  ...  ,  and  n  =  0,  1,  ...  ,  the  matrices 
^/),  R'„{t)  and  ejit)  can  be  written 
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(2.23) 


KU)  0  1 
0  ;//J  ■ 

mxm 

Also,  the  initializations  in  (2.24)  and  (2.25)  are  straightforward  from  (2.i2)-(2.21).  The  residual-error  lattice 
in  Section  III.B  of  [3]  then  becomes  Algorithm  2.1  for  Problem  2.1. 

Initialization  2.1  (For  Unwindowed  AR  Problem) 


pxp 

II 

o 

_E 

1 _ 

II 

o 

[  0  /  ]• 

(2.24) 

mxm 

A:„(0  =  0  ,  n  >  t  ^  0 . 

(2.25) 

Algorithm  2.1  (Residual  Error  Lattice  for  the  Unwindowed  .AR  Problem) 
For  each  /  >  1, 


^o(0  =  [eo(0  I  0]  =  CKO  I  0],  (2.26) 

^'(0  =  .y  =  ro^(0?o(0  +  1),  Go(0  =  /.  (2.27) 

For  n  =  0  to  min{t  -  (..S')  (where  N  is  the  ma-iamum  order), 

=  ^-A;^i('-  1)  +  ^I(i)G;{t-  l);„(r-  1)  (2.28) 

G.+  i(0  =  G„it)  -  ?„(i)  ^;,%)?J'(t)  (2.29) 

<+,(0  =  Ko-  1)  -  /:„+,(0C(0/f„+,(0  (2.30) 

A„%,(0  =  ^^0  -  /f^,(o^„^(^-  i)a:„+i(0  (2.31) 

en+,(0  =  e^t)  -  ?„(t-  I)A:J^,(0  (2.32) 

^+i(0  =  r„(l-  1)  -  ^K0/?f(0/:,+,(0.  (2.33) 

As  in  [3],  Gj(i)  is  an  m  x  m  matrix.  For  a  matrix  M,  M*  means  any  matrix  such  that  MM*M  =  M. 


The  residual-error  lattice  for  the  unwindowed  AR  problem  has  the  same  form  as  the  residual-error  lattice 
in  [3]  for  the  prewindowed  problem,  but  in  the  lattice  for  the  unwindowed  problem  the  anays  KJ^t),  r^i) 


enU)  =  [e.(0  I  0  ]  ■  ^nU)  = 


mxp  mxm 


m 

0 


pxp 
mx  p 


1 .  ht)  =  [ 
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and  R'„{i)  are  larger  than  the  corresponding  arrays  in  the  prewindowed  lattice.  WTrile  the  matrices  /?'(/)  and 
G„{t)  have  the  same  respective  dimensions  as  in  the  prewindowed  problem,  the  matri.x  R^{i)  in  the  unwin- 
dowed  problem  has  dimension  (p  +  m)  x  (p  +  m)  as  opposed  to  p  x  p  for  the  corresponding  matrix  in  the 
prewindowed  problem.  The  unwindowed  lattice  should  be  more  complex  than  the  prewindowed  lattice, 
but  In\erting  the  larger  matrix  R'j^t)  is  both  undesirable  and  unnecessary. 

From  the  way  in  which  R^(t),  R^t)  and  KJit)  arise  in  [3]  (they  arc  used  in  computing  projections),  it 
follows  that  ^(R!Xt))  and  ,^(A!^,(f))c:  1)).  According  to  Lemma  A.l  in  the  Appendix 

then,  we  can  generate  an  R^{i)  directly  with 

5:1(0  =  1)  +  51/ -  i)/^„';,(oC,(o/^„+,(o5l/- 1) .  (2  30') 

5io)  =  =  5(0) . 

where  M'  is  the  usual  pseudo  inverse  of  a  matrix  .V/  [9].  For  the  initialization  in  (2.30'),  recall  (2.24).  The 
y?,‘(f)  generated  by  (2.30'i  is  not  /?r(0  in  general. 

In  the  most  efficient  version  of  Algorithm  2.1  then,  (2.30')  replaces  (2.30)  and  R^{i)  is  used  for  R^{t)  in 
(2.30).  (This  is  true  for  sohing  the  FIR  problems  in  Sections  3  and  4,  also.)  Although  it  is  not  necessary  , 
it  is  most  natural  to  use  G'.{t)  in  (2.28). 

.Also,  we  normally  do  not  compute  /C*(0  and  G~(t)  exactly  when  these  matrices  are  singular  or  near 
singular.  Rather,  we  choose  a  small  number  S  and  for  det(/?'(/))  <  5  or  det(G„(0)  <  <5,  we  use  the  approxi¬ 
mation 

(x/  -f-  (2.34) 

for  small  positive  x.  Numerical  studies  in  [10]  show  that  using  (2.34)  to  approximate  /C*(0  and  G*(0 
produces  negligible  error  in  the  lattice  results.  For  the  numerical  results  in  Section  5,  we  used 
a  =  3=  lO-'o  . 

The  AR  coefficients  /!.,,  ,(/)  for  problem  2.1  can  be  generated  at  any  i  by  the  following  algorithm.  Like 
Algorithm  2.1,  Algorithm  2.2  is  derived  from  the  corresponding  algorithm  in  [3]  by  embedding  Problem 
2. 1  in  Problem  2.2  and  then  eliminating  the  parts  of  the  expanded  prewindowed  algorithm  that  are  not  need 
for  the  solution  of  Problem  2.1. 
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Algorithm  2.2  (AR  Coefficients  for  the  L’nwindowed  AR  Problem) 
For  n=  I,  2,  l  —  I  and y  =  1 ,  2,  ...  ,  n 


Q+i  =  C,  j(t)  -  .  (2.35) 

LSnjiO-  C„Jt)GUt)7„(i)]  -  A„Jt)R^(i)K„^,(i),  (2.36) 

'<n+Uj(l)  =  ^njiO  -  LB„jin  -  C„  Jc)G^(t)7„((n  K%-  l)^n+,(f)  .  (2.37) 

The  end  conditions,  for  n  =  0,  1,  ...  ,  /  —  1,  are 

/f.+  ,^+,(f)  =  IKV-  l)\^~  (2.38) 

=  R^it)K„a‘),  (2.39) 

C„^,^^,(f)=  [^nt)]^,~;j(/).  (2.40) 

where  means  the  top pxp  block  in  ^{t). 
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3.  The  TMR  Problem  (Joint  IVocess  Estimation) 

Now  we  assume  that,  in  addition  to  the  sequence  oi  my.  p  measurement  matrices  ;-(0.  ba\e  a  se¬ 
quence  of  m  ■  <7  m.easurement  matrices 


=  [x'(0-v\0  ...  x-'Cr)]. 


By  an  n'*-order  FIR  model,  we  mean 


^0  =  +  e„(0. 

J=i 


(3.1) 


(3.2) 


where  the  A„  .'s  are  /»x  <7  matrices  referred  to  as  FIR  coefficients,  and  the  /  *  column  of  A^j  is  denoted  by 

A 

Ai^j.  We  denote  the  history  of  x(i)  by 


0(0  =  [o'(0d>‘(0  ...  =  [^^(Ox^(t-  i)x^(f-2)  ...  f  . 


(3.3) 


Problem  3.1.  (Unwindowed  FIR  Problem)  For  /  =  1,  2,  ...,  and  n  =  0,  1,  ...  .  r,  compute 
P,,,.,-,(^+  l)<?(0  =  l)0,-i-n<^(O;  equivalently,  (for  n=  1,  2,  ...  ,  i)  find  matrices  .■5,./0  such 

that,  for  each  i  (/  =  1.  ...  ,  <7),  the  /■'*  columns  of  these  matrices  minimize 


1  x\x)  -  ^><r  -y+  l)^'.,.(0  1' 


(3.4) 


M 


Problem  3.2.  (Prewmdowed  FIR  Problem)  Assume  x(t)=  0  and  y{l)-=  0  for  f  <  0.  For  f  =  0,  1,  2,  ..., 

and  n  =  0,  1 . t+  1,  compute  P^l  +  1)0(0;  equivalently,  (for  n=  1,  2,  ...  ,  ^ -I-  1)  find  matrices 

An,j(0  such  that,  for  each  /(/  =  1,  ...  ,  <7),  the  /'*  columns  of  these  matrices  minimize 


J^l)  =  I  JcV)  -  -J+  l)^^.y(0l' 


(3.5) 


x*0 


J-1 


As  in  Section  2,  we  embed  the  unwindowed  problem  in  the  prewindowed  problem  by  defining  the 
measurement  sequence  y(t)  for  the  prewindowed  problem  in  terms  of  the  measurement  sequence  ><0  for  the 
unwindowed  problem  according  to  (2.19)-(2.21).  We  have 


1)0(0  =  ^Pn,.,^-n0+  1)0(0  0].  /  =  I.  2 .  n  =  1,  2 . 


(3.6) 


9 


so  that  the  first  (/  +  1  —  n)m  rows  of  the  first  q  columns  of  PJ^t  +  1)<^(0  arc  equal  to  the  first  (t  +  1  —  n)m 
rows  (the  only  possibly  nonzero  rows)  of  +  1)<^(/).  Therefore,  the  solution  to  Problem  3.1  for  the 

sequences  x(t)  and  is  embedded  in  the  solution  to  Problem  3.2  for  the  sequences  x(t)  and  ji(/) . 

To  obtain  a  lattice  for  Problem  3.2,  we  must  append  some  new  ujxiatc  equations  to  the  prewindowed 


lattice  in  [3].  For  nonnegative  integers  /  and  n,  we  define 

fnio  =  irM  Z'^(o]  =  [/  -  hi  + 1)]  m .  ot) 

A  ^ 

e„{t)— top  m  rows  of  f„{t),  (3.8) 

^n+i(0=  matrix  whose  (i,  y)  element  is  <f„it),  6^(f)>  ,  (3.9) 

A  A  .  A  . 

R„it}=  q^  q  matrix  whose  (i,  j)  element  is  7^(0)  .  (3.10) 


where,  as  in  Section  2,  P„{t)  is  the  orthogonal  projection  onto  H„(t)  when,  in  (2.7),  p  is  replaced  by  p  and 
i>\t—j)  is  replaced  by  ip‘{t—  j)  for  i  =  1,  ...  ,p. 

To  derive  order  updates,  we  define 

P^(t)  =  orthogonal  projection  onto  span{b),{t)  ...  b^(t)}  ,  (3.11) 

with  b'jit)  defined  as  in  (2.13).  We  have  then 

[/-P.+  ,(0]=[/-P.V-  I)]Cf- WJ.  (3.12) 

From  (3.7),  (3.9),  (3.12)  and  (2.18),  we  obtain 

L,io  =  frO)  -  Pnimo  =  m  -  ■  (3.13) 

It  is  straightforward  to  use  (3.7)-(3.10)  and  (3.13)  to  derive  Initialization  3.1  and  the  order  updates  in  Al* 

A 

gorithm  3. 1 .  With  the  foregoing  equations  in  this  section,  the  derivations  of  the  time  update  for  KJt)  and 
of  Algorithm  3.2  are  similar  to  those  in  Sections  III.B  and  III.C  of  [3]  (sec  [10]).  The  most  difficult  of 

A 

these  derivations  is  that  for  the  time  update  of  Kft). 
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Initialization  3.1  (For  the  Unwindowed  FIR  Problem,  append  this  initialization  to  Initialization  2.1.) 


^(0)  =  x^(0)x(0).  ^,(0)  =  [  0  1  x^(0)] 

(3.14) 

^n+i(0  =0-  n  >  '  ^  0 

(3.15) 

Algorithm  3.1.  (For  the  Unwindowed  FIR  Problem,  append  this  algorithm  to  Algorithm  2.1.) 

For  t  >  1 

e^t)  =  x{t) ,  ^'(0  =  x\t)jit)  +  A  ^'(/-  1) .  (3.16) 

For  n  =  0  to  / 


=  ;.A;^,(:-  1)  +  eJ'(t)G;^(t)?„(t) 

(3.17) 

(3.18) 

=  k(0  -  -^„+,(f)C(0^„+,(0. 

(3.19) 

Algorithm  3.2.  (For  the  FIR  problem,  append  this  algorithm  to  Algorithm  2.2.) 

For  n  =  1,  2,  ...  ,  /  and  y  =  1,  2,  ...  ,  n. 

(3.20) 

The  end  condition  .  for  /i  =  0,  1,  ...  ,  r,  is 

(3.21) 

Algorithms  3.1  and  3.2  with  Initialization  3.1  solve  Problem  3.1.  When  the  solution  to  Problem  3.2 
with  arbitrary  ^/),  /  ^  0  ,  is  desired.  Algorithm  3. 1  can  be  appended  to  the  prewindowed  lattice  in  [3].  For 

A  ^ 

this  case,  the  only  change  in  Initialization  3.1  is  that  1^,(0)  =  x^0)y(0)  must  be  used  in  (3.14).  Algorithm 
3.2  generates  the  full  matrices  A(t)  for  Problem  3.2  if  A^j(t\  B^/t)aDd  [/?,  (0],.;  are  replaced  by 
/Oand  ^(0.  where  is  the  matrix  generated  by  the  algorithm  in  [3]  for  the  AR  coefficients. 
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4.  Including  Initial  Parameter  Estimates 


The  algorithms  and  initializations  in  the  previous  sections  do  not  use  initial  parameter  estimates.  To 
include  initial  estimates  of  the  FIR  parameters  corresponding  to  a  given  order  N,  we  need  only  change  In¬ 
itializations  2.1  and  3.1. 

Problem  4.1.  (Unwindowed  FIR  Problem  with  Initial  Parameter  Estimates)  For  r  =  1,  2,  ...,  and 

A 

n=  1,  ...  ,  A’,  find  matrices  such  that,  for  each  /(/  =  1 . ^),  the  /“*  columns  of  these  matrices 

minimize 

n 

-I-  1  Z.jit)  -  A^s.j  Id  .  (41) 

where  Aff  j  is  the  initial  estimate  of  A‘^j  ,  D  =  diag  {5,,  5j.  ...  ,  5,}  ^  0, 1  a  |^  =  a^Da  and  n  >  0. 

Like  Problems  2.1  and  3.1,  Problem  4.1  can  be  embedded  in  a  prewindowed  problem.  We  defme  a 
prewindowed  FIR  problem  similar  to  Problem  3.2  except  that  the  nonzero  data  begins  at  r  =  —  p,\  instead 
of  :=  0.  For  j=  2,  ...  ,  .V  and  k=  I,  2,  ...  ,  p  ,  we  define 


row  k  of  A 

0  0  ...  0 


k*  place 


*9 

\{m-\)xq 


1/2  cl/2 


00 


0 


00 

00 


llxp 

J{/n-l)x 


(4.2) 


1  -  Xk)  =  a, v.yjc.  I  :<  y  ^  A',  I  <  k  ^  p. 


(4.3) 


U'e  replace  (2. 19)  with 


An  = 


[O  . 


t  =  —kN ,  \  <  k  <  p , 

all  other  f  <  0. 


(4.4) 


For  t^  0,  y{t)  is  defmed  by  (2.20)-(2.21).  It  is  straightforward  to  see  that  the  prewindowed  FIR  problem 
for  these  histories  of  jc(r)  and  is  equivalent  to  Problem  4.1.  Furthermore,  for  t  =  0,  the  various  matrices 
in  (2.12)-(2.18)  and  (3.8)-(3.10)  can  be  computed  directly  from  the  definitions  of  x{t)  and  AO  yield  the 
following. 
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Initialization  4. 1 


roiO)  =  [  0  I  /  ] .  ;„(0)  =0  n=l.2,  ,  N-  1  (4.5) 

pxp  pxp 

"^(0)=  jj,  <(0)=  gj,  n=\,2 . V-l  (4.6) 

mx/n  mxm 

^„+i(0)=0,  n=0.  1 . N-2  (A.l) 

G^O)  =  I,  G„(0)  =0.  n=l,2 . N-  1  (4.8) 

,v 

^o(0)  *  X  ^  (4-9) 

./=! 


^,(0)  =  I  x^(0)],  .^„^,(0)  =  I  0],  «=  i,  2,  ...  ,  a  -  i  .  (4.10) 

For  the  solution  to  Problem  4.1,  Initialization  4.1  replaces  Initializations  2.1  and  3.1.  Algorithms  2.1, 
2.2,  3. 1  and  3.2  are  used  as  they  are  stated,  except  that  the  maximum  n  in  Algorithms  2. 1  and  2.2  is  .V  —  2 
and  the  maximurn  n  in  Algorithms  3.1  and  3.2  is  N  —  1. 

Note  that  for  =  0,  Initialization  4. 1  reduces  to  Initializations  2. 1  and  3. 1 .  On  the  other  hand,  for 
Affj=  0,  Initialization  4.1  is  different  from  the  prcwus  initializations  because  they  imply  no  initial  param¬ 
eter  estimates,  which  is  different  from  setting  initial  parameter  estimates  equal  to  zero. 
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5.  Example 


We  used  a  sixth-order  AR  model  with  m  =  p  =  1  to  generate  the  data  sequence  }it),  /  =  1,  2,  ....  100, 
from  nonzero  initial  conditions.  The  true  parameters  and  initial  conditions  are  given  in  Tabic  1.  These 
AR  coefficients  represent  a  single  measurement  from  a  system  of  coupled,  lightly  damped  oscillators  sam¬ 
pled  at  the  rate  of  50  Hz. 

TABLE  1 

True  AR  Coefficients  and  Initial  Data 


J 

Ml -7) 

1 

4.0021 

0.27 

2 

-7.8919 

0.18 

3 

9.6013 

0.029 

4 

-7.5673 

0.096 

5 

3.6927 

0.35 

6 

-0.89871 

0.50 

First,  we  will  compare  results  obtained  with  the  unwindowed  AR  lattice  in  this  paper  and  results  ob¬ 
tained  with  the  prewindowed  lattice  in  [3].  Rather  than  show  the  estimates  of  the  AR  parameters,  we  plot 
in  Figure  1  the  three  frequencies  computed  from  the  parameter  estimates  according  to 

cu  =  Im{{  logC)x  50) ,  (5.1) 

w  here  (  is  an  eigenvalue  of  the  AR  model.  The  frequencies  corresponding  to  the  true  parameters  are 

cuj  =  14.549,  cuj  =  40.743,  cuj  =  58.862.  (5.2) 

Also,  we  plot  in  Figure  2  the  one-step-ahead  prediction  error;  at  any  time  r,  this  is  the  error  between  y{i) 
and  the  summation  on  the  right  side  of  (2.2)  evaluated  with  the  estimated  parameters. 

In  Figures  1  and  2,  solid  curves  correspond  to  parameters  generated  by  the  algorithms  in  Section  2  of 
this  paper,  and  the  dashed  curves  correspond  to  p>arameteTS  generated  with  the  prewindowed  lattice  in  [3]. 
In  both  lattices,  we  used  the  forgetting  factor  A  =  .98.  All  numerical  results  presented  here  are  for  n  =  6. 
If  enough  data  is  used,  the  frequencies  corresponding  to  the  parameter  estimates  from  the  prewindowed 
lattice  eventually  will  converge  to  the  true  frequencies.  This  convergence  is  faster  for  smaller  values  of  A, 


but  to  make  the  convergence  significantly  faster,  >.  must  be  so  small  that  the  lattice  filter  becomes  very 
sensitive  to  noise. 

To  illustrate  the  effect  of  setting  initial  parameter  estimates  with  Initialization  4.1.  we  added  zero-mean 
white  noise  to  the  data  generated  by  the  AR  model  for  30<  /  <  50.  The  standard  deviation  of  the  noise 
was  .01,  which  is  between  1.5“  'o  and  2'’'o  of  the  amplitude  of>fr).  We  generated  two  sequences  of  parameter 
estimates  with  the  unwindowed  FIR  lattice  (i.e..  Algorithms  2.1,  2.2,  3.1  and  3.2)  and  Initialization  4.1  with 
jV  =  6.  We  obtained  an  FIR  problem  by  setting  jc(/)  =  •+  1),  /  ^  1  •  For  the  first  sequence  of  parameter 

estimates,  we  started  at  r  =  1  with  /i  =  0  (i.e.,  no  initial  parameter  estimates)  and  continued  until  r=  100. 
For  the  second  sequence  of  parameter  estimates,  we  also  started  at  /  =  1  with  /i  =  0,  but  at  r  =  30,  w  e  re¬ 
initialized  the  FIR  lattice,  using  ^  =  100  ,  Z)  =  /  and  Af,  j  =  As  j(29)  in  Initialization  4.1. 

For  each  of  these  sequences  of  parameter  estimates.  Figure  3  shows  the  norm  of  the  error  between  the 
estimated  parameter  \  ector  and  the  true  parameter  vector.  For  1  <  /  S  29,  the  two  sequences  of  parameter 
estimates  are  identical. 

For  l  >  30,  the  measurement  noise  affects  the  parameter  estimates  generated  with  no  resetting  much 
more  than  it  affects  the  parameter  estimates  generated  after  resetting,  since  the  reinitialized  FIR  coefficients 
are  held  near  the  correct  values  at  t=  29  by  the  heavy  weighting  /i  .  With  smaller  forgetting  factor  the 
parameter  estimates  would  return  to  true  values  faster  after  i  =  50,  but  smaller  /.  would  defeat  the  purpose 
of  n  . 

We  compared  the  performance  of  the  lattices  in  this  paper,  the  prewindowed  lattice  in  [3],  and  the  un¬ 
windowed  lattice  in  [7]  on  the  current  example.  Since  [7]  did  not  indicate  how  to  generate  estimates  of 
either  AR  or  FIR  parameters,  we  compared  the  values  produced  by  our  lattices  for  the  forward  residual 
errors  e^t)  and  ejii)  to  the  corresponding  quantities  produced  by  the  lattice  in  [7].  We  made  this  compar¬ 
ison  for  the  AR  problem  correspondmg  to  Figures  I  and  2,  so  that  e^/)=  e«(/+  1).  The  quantities  ejii), 
ejii)  and  /?'(f)  should  correspond,  respectively,  to  the  quantities  e^,  and  (with  A:  =  f)  in  [7]. 

The  Exact  Initialization  in  [7],  with  the  nonzero  value  for  jc,  indicated  in  [7],  >ielded  values  for  the 
forward  error  ,  the  error  and  squared  error  norm  all  very  close  to  the  corresponding  quan¬ 
tities  from  the  prewindowed  lattice.  Both  the  lattice  in  [7]  and  the  prewindowed  lattice  yielded  values  of 

Sind  on  the  order  of  10"’  and  values  of  on  the  order  of  !()■'.  Our  unwindowed  lattice  yielded 
values  of  e/^t)  and  e^r)  on  the  order  of  KT's  and  values  of  R^t)  on  the  order  of  1(T'*.  These  results  indicate 
that  our  unwindowed  lattice  was  soKing  the  unwindowed  problem  and  that  the  lattice  in  [7]  was  not. 
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When  we  shifted  all  data  for  A  >  0  for^v  ard  by  one  step  and  set  jCj  =  0  in  the  Exact  Initialization,  the 
lattice  in  [7]  yielded  and  that  agree  with  the  corresponding  quantities  from  the  prewin¬ 

dowed  lattice  to  nine  digits,  all  that  we  printed.  We  also  ran  the  lattice  in  [7]  with  the  Soft-Constraint  In¬ 
itialization  on  page  368  in  [7]  and  compared  the  results  to  those  from  our  lattice  with  soft-contraint 
initialization.  .Again,  the  for^vard  errors  from  the  lattice  in  [7]  were  much  larger  than  those  from  our  lattice. 

Our  analysis  of  the  derivations  and  algorithms  in  [7]  indicates  that,  when  all  data  is  zero  for  r  <  0,  the 
lattice  in  [7]  with  the  Exact  Initialization  solves  the  prewindowed  problem,  but  the  lattice  in  [7]  with  either 
initialization  in  [7]  does  not  solve  the  un windowed  problem.  These  conclusions  are  confumed  by  our 
numerical  results. 

6.  Conclusions 

The  purpose  of  the  paper  has  been  to  derive  efficient  lattice  filters  to  solve  Problems  2.1,  3.1  and  4.1. 
Problem  2.1,  the  unw  indow  ed  .AR  problem  with  no  initial  parameter  estimates,  is  solved  by  Algorithms  2.1 
and  2.2  with  Initialization  2.1.  Tlie  residual-error  lattice.  Algorithm  2.1,  is  recursive  in  both  time  and  order. 
Algorithm  2.2.  which  generates  the  estimates  of  the  AR  coefficients,  also  has  a  lattice  structure,  but  this 
algorithm  is  recursive  in  order  only.  It  can  be  run  at  any  time  t,  but  it  need  not  be  run  at  every  t.  For 
Problem  3.1,  the  unwindowed  FIR  problem  with  no  initial  parameter  estimates.  Algorithms  3.1  and  3.2  and 
Imtialization  3.1  are  run  along  with  Algorithms  2.1  and  2.2  and  Initialization  2.1.  The  solution  to  Problem 
4.1,  the  unwindowed  FIR  problem  with  initial  parameter  estimates,  also  requires  Algorithms  2.1,  2.2,  3.1 
and  3.2,  but  Irutiahzation  4  1  replaces  Initializations  2.1  and  3.1. 

To  make  the  lattices  in  this  paper  most  efficient,  the  generalized  inverse  R^(t)  is  generated  directly  with 
(2.30  );  RJt)  need  not  be  generated  and  (2  30)  need  not  be  used.  Lemma  A  l,  on  which  (2.30')  is  based,  is 
a  generahzation  of  a  well  known  formula,  whose  diverse  appheations  were  surveyed  recently  in  [11]. 
Lemma  A.  1  appears  to  be  new. 
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Figure  1 ;  Frequency  Estimates 


Figure  2;  One-Step-Ahead  Predition  Erron 
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log  (Norm 


Figure  3:  Norms  of  Parameter  Errors 
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APPENDIX 


Lemma  A. I  Let  A,  B  and  C  be  real  p  /p.  py-p  and  py.p  matrices,  respectively,  with  ^iB)<z  ^i(A)  , 
^(Q.  .V(/l)c  \(B^}  and  .V(C)c  .V(S).  Suppose  that  matrices  A .  C.  ,  C.  .T’  and  C  satisfy 


A  =  A  -  BC""B^.  (.-I.l) 

"C^C-B^A'^B.  (A.2) 

A""  =  A""  +  A""  BC^  B^A"" ,  (.-1.3) 

and 

AA  ^A  =  A,  CC^C  =  C .  CC^^C  =  C.  (.-1.4) 

Then 

A  A"  A  =  A  .  (.-1.5) 


Proof  The  inclusion  \(A)c:  .\(B^)  and  .-l[.-l'.4  —  /]  =  0  imply 

B'^A  A  =  B^.  (.-1.6) 

Then 

B ^A  "A  =  B ^A  =.-1  -  B '’’a  ""BC^B ^  =  CC^B ^-B'^'a  ^BC"B ^  =  CC"B ^ .  (.4.7) 

Hence 

^(5^.4=7)c,5?(C).  (.4.8) 

.Also,  replacing  A,  A',  B,  C  and  (f  in  (.A. 7)  and  (A.8)  with  C.  Cf,  B^ .  A  and  A' ,  res’^ctively,  yields 

BC^C^AA^B.  (.4.9) 

Now, 


/1[.4“- .4''5C''fl^.4=].4  =AA'‘A  -f  ( .4 /I “5)C^fi ^.4 N4 

=  .7.4  =.T  +5C=CC^5^.4'’7  =  .7.4  "'.7  +  BC=fl74N7  =  .4.4  n7  =  .7  . 


The  last  identity  in  (A.lOj  follows  from  i?(.4  )c  ^(A) ,  which  follows  from  ^(B)c.  ^{A)  and  (A.li. 


In  general.  .4’  ^  A  (the  pseudo  inverse  of  .4  ),  even  when  .4  and  C  are  used  for  .4’  and  C,  respecti\ely, 
in  the  right  side  of  (A. 3).  For  example,  if 


then 


and 


o]-  C^CA^O. 

A"  =  A~'  +  A~'B^B'^A~'  =  A~'  =  [  Jj  ~ ‘  ] 


(.4.12) 


(.-1.13) 
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ABSTRACT 

A  digital  adaptive  controller  for  a  robotic  manipulator  with  a  sliding 
flexible  link  is  presented.  The  most  important  feature  of  the  controller  is  its 
capability  to  vary  the  order  of  the  control  law  adaptively.  This  capability 
results  from  using  a  lattice  filter  for  adaptive  parameter  estimation.  The 
superiority  of  the  variable-order  adaptive  controller  to  a  fixed-order  adaptive 
controller  Is  demonstrated  by  numerical  simulations.  In  which  the  manipulator  Is 
represented  by  the  nonlinear  equations  of  motion  for  a  finite  element  model  of 
the  manipulator. 
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1.  INTRODUCTION 

Adaptive  control  of  robotic  manipulators  with  rigid  links  has  been  studied 
by  many  authors  [BDSl,  DDl,  K2,  KGl.  LEI,  NTl,  THl,  VKl].  While  several  authors 
have  studied  nonadaptlve,  usually  optimal,  control  of  manipulators  with  flexible 
links  [BMWl,  CSI,  F2,  UNMl],  the  literature  contains  only  limited  treatment  of 
adaptive  control  of  manipulators  with  flexible  links.  In  [NMl,  NMBl],  an 
unknown  payload  was  estimated  on-line  to  update  optimal  control  gains  used  to 
control  a  linear  model  of  a  flexible  link.  Parameters  were  Identified  in  [RCl] 
and  then  used  to  design  a  steady-state  Ilnear-quadratIc-Gausslan  compensator, 
which  was  used  to  control  a  flexible  one-link  manipulator.  In  [CLl],  a  static 
elastic  deflection  was  modeled  In  one  simulation  of  adaptive  pole  placement  for 
the  Stanford  arm.  A  discrete-time  adaptive  controller  designed  for  a  rigid  link 
was  applied  to  a  flexible  link  In  [Yl],  and  the  adaptive  controller  did  not 
appear  to  suppress  all  oscillations  of  the  link  about  the  final  position. 

This  paper  presents  a  digital  adaptive  controller  for  a  manipulator  with 
two  links,  the  second  of  which  Is  a  flexible  beam  that  slides  out  of  the  first 
link  (I.e.,  there  Is  a  prismatic  joint).  The  simulation  model  of  the  manipula¬ 
tor  Is  nonlinear  with  time-varying  coefficients,  but  the  adaptive  controller  Is 
based  on  a  linear  Auto  Regressive-Moving  Average  (ARMA)  model  of  the  Input/out¬ 
put  characteristics  of  the  plant.  Both  the  parameters  and  the  order  of  this 
ARMA  model  vary  adaptively.  Becasue  the  control  law  is  based  explicitly  on  the 
ARMA  model,  the  adaptive  control  algorithm  is  Indirect. 

The  most  Important  feature  of  the  adaptive  controller  In  this  paper  is  that 
Its  order  can  vary  adaptively.  This  capability  is  Important  in  the  problem  here 
because  the  flexible  link  abruptly  ceases  sliding  and  the  associated  large  axial 
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deceleration  has  the  effect  of  a  lateral  Impulse  on  the  link.  The  Impulse  ex¬ 
cites  previously  unexcited  elastic  modes  of  vibration,  so  that  the  effective 
order  of  the  plant  Increases.  Such  changes  In  plant  order,  which  might  result 
also  from  Impacts  or  releasing  payloads,  are  handled  much  better  by  a  controller 
whose  order  can  vary  adaptively. 

The  order  of  the  controller  can  vary  adaptively  because  the  parameter  esti¬ 
mator  Is  a  least-squares  lattice  filter,  which  Is  an  algorithm  for  least-squares 
parameter  estimation  that  Is  recursive  In  both  time  and  order.  The  lattice 
filter  Is  more  efficient  than  the  standard  least-squares  method  for  large 
orders,  and  It  Is  numerically  stable  [LLl].  Lattice  filters  have  become  promi¬ 
nent  In  adaptive  signal  processing  [Fl,  GSl,  HMl,  LFMl,  LMFl,  LSI].  Their  use 
In  control  and  Identification  of  flexible  structures  has  been  studied  In  [MSI, 
SMI,  SM2,  SM3,  Wl,  WGl],  and  recently  the  capability  of  high-order  lattices  to 
Identify  many  modes  of  very  flexible  structures  has  been  demonstrated  In  [Jl, 
JGl,  JG2]. 

Recursive  (In  time)  least-squares  estimation  Is  used  widely  In  adaptive 
parameter  Identification.  See  [LSI,  GSl]  and  many  other  references  for  the 
classical  algorithm  and  Its  convergence  properties.  This  method  has  one  serious 
limitation  for  Identification  of  systems  like  the  manipulator  in  this  paper:  the 
classical  recursive  least-squares  algorithm  Is  based  on  a  fixed-order  input/ 
output  model.  In  a  flexible  system,  different  numbers  of  inodes  may  be  excited 
at  different  times,  especially  In  problems  with  Impulse  phenomena.  For  such 
problems,  determination  of  effective  plant  order  Is  needed  along  with  iden¬ 
tification  of  parameters.  With  Its  order-recursive  property,  the  lattice  filter 
can  Identify  the  number  of  substantially  excited  inodes  of  a  flexible  system  as 
well  as  the  parameters  of  a  digital  input/output  model. 
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Section  2  describes  the  manipulator,  the  finite  element  model  used  for 
simulations,  and  the  ARMA  model  on  which  the  control  law  Is  based.  Section  3 
discusses  the  adaptive  parameter  estimation  and  presents  the  lattice  filter 
algorithm.  The  adaptive  control  law  Is  discussed  In  Section  4,  along  with  the 
reference  signal,  learning  period,  control-torque  bounds  and  the  criterion  for 
changing  the  order  of  the  controller.  Numerical  simulations  are  presented  and 
discussed  In  Section  5.  and  the  paper's  conclusions  are  stated  In  Section  6. 

2.  MANIPULATOR  MODEL 

Figure  2.1  shows  the  manipulator  to  be  controlled,  which  consists  of  the 
rotor  whose  center  Is  fixed  at  point  0,  the  rigid  link  Mj  and  the  flexible  link 
^2^  The  rotor  Is  modeled  as  a  rigid  disc  and  the  rigid  link  Is  cantilevered  to 
the  rotor.  The  flexible  link  Is  a  uniform  Euler-Bernoulll  beam  that  slides  in 
the  rigid  link.  The  free  end  of  the  flexible  link  carries  a  payload,  modeled  as 
a  point  mass  Mp|^.  As  Table  2.1  Indicates,  the  manipulator  Is  unusually  long  and 
slender.  This  design  Is  motivated  by  possible  space  applications  of  robotics 
and  by  the  desire  for  a  highly  flexible,  rather  unwieldy  manipulator  to 
challenge  the  controller. 

We  assume  that  the  radial  motion  of  the  flexible  link  Is  controlled  by  an 
actuator  with  sufficiently  wide  bandwidth  and  short  time  constant  to  make  the 
flexible  link  follow  a  specified  radial  position  profile  r(t)  exactly. 
Consequently,  we  treat  r(t)  and  Its  derivatives  as  time-varying  parameters.  The 
control  variable  In  this  paper  Is  a  torque  u  that  acts  on  the  rotor  at  0. 
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TABLE  2.1 


Manipulator  Data 


Combined  Rotor  and  Rigid  Link 

Radius  of  rotor  plus  length  of  rigid  link  ^  L  ^  3m 
Combined  moment  of  inertia  about  rotor  axis  >  168.48 

Rexible  Link(steel) 

Length  ^  3m 

Length  of  segoent  outside  of  rigid  link  >  r  (Im^r^lm) 
Cross  section  •  0.02  mx 0.02  m 
Mass  per  unit  volume  ■  7.8x10^ 

Modulus  of  elasticity  ■  E  *  200x10*  Nfm^ 

Voigt-Kelvin  damping  coefficient  *  *  0.001 

Fundamental  cantilevered  beam  frequency* 

For  r—  1  m,  /j  —  16.336 Hz 
For  r—2/n,  /,  —  4.089 //z 


In  the  dynamic  model  of  the  manipulator: 

•  first-order  (linear)  transverse  deformations  of  the  flexible  link  are 
modeled; 

•  axial  elastic  deformations  are  neglected; 

•  coupling  between  rigid-body  angular  velocity  and  elastic  displacements  is 
Included; 

•  the  contribution  of  the  Inertial  axial  load  to  bending  stiffness  Is 
Included; 

•  torques  due  to  gravity  are  Included. 

For  simulating  the  response  of  the  manipulator,  we  use  three  finite  ele¬ 
ments  of  equal  length  to  approximate  the  part  of  the  flexible  link  outside  the 
rigid  link.  Since  this  portion  of  the  flexible  link  varies  with  time,  the 
length  r(t)/3  of  each  of  the  three  elements  must  vary  with  time.  We  use  cubic 
B-splInes  [SI]  as  basis  functions.  This  means  that  we  have  three  elastic 
degrees  of  freedom,  which  we  take  to  be  the  transverse  elastic  displacements  of 
nodes  2,  3  and  4.  (Node  1  Is  the  point  on  the  flexible  link  at  the  end  of  the 
rigid  link;  node  4  Is  the  end  of  the  flexible  link  to  which  the  payload  Is 
attached.)  In  all  then,  there  are  four  degrees  of  freedom  in  our  simulation 
model  of  the  manipulator.  We  represent  the  rigid-body  degree  of  freedom  by  the 
angle  6. 

For  the  finite  element  model  of  the  manipulator,  the  generalized  displace¬ 
ment  vector  is  q  ■  [6  q2  (>3  where  q2>  q3  and  q^  are  the  transverse  elastic 
displacements  of  nodes  2-4  on  the  flexible  link.  Lagrange's  equations  for  the 
finite  element  model  have  the  form 

M(t)ii  +  D(t)q  +  K(t)q  ■  Qg(t)g  sin  6  +  Bu,  (2.1) 
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where  the  4x4  matrices  M(t),  D(t)  and  ic(t)  and  the  4-vector  Og(t)  have  the  fonr.s 


M(t)  -  M(q.r(t)).  (2.2) 
D(t)  -  Dj(q.q.r(t))  +  r(t)02(q.r(t)) .  (2.3) 
>C(t)  -  Ki(q.q.r(t).g)  +  +  r(t)IC3(r(t)).  (2.4) 
Og(t)  -  Og(r{t)),  (2.5) 
B  -  [1  0  0  0]*^.  (2.6) 


and  g  Is  the  acceleration  of  gravity.  Of  course,  the  generalized  coordinates  q 
vary  with  time,  but  we  want  to  emphasize  that  r(t)  Introduces  time-varying  coef¬ 
ficients  Into  the  equations  of  motion.  We  model  small  Voigt-Kelvin  viscoelastic 
damping  (See  [CPI])  In  the  link,  which  means  that  the  matrix  D(t)  contains 
CqKjj  where  Is  the  part  of  K(t)  that  represents  the  bending  stiffness  of  the 
flexible  link  and  Cq  Is  a  damping  coefficient.  For  a  complete  derivation  of  the 
equations  of  motion,  see  [Kl]. 

The  flexible  link  slides  out  of  the  base  link  according  to  the  length- 
versus-tlme  profile  In  Figure  2.2.  During  the  time  between  0  and  tQ,  which  will 
be  used  for  preliminary  parameter  Identification  by  the  adaptive  controller,  the 
radial  velocity  r  Is  zero  and  the  rigid-body  angle  remains  near  zero  (See 
Section  4.3).  Between  tg  and  t^,  the  flexible  link  has  the  circular  radial 
motion  profile  with  positive  f  and  f  indicated  In  Figure  2.2.  Between  and 
t2,  f  Is  constant,  and  starting  at  t2,  r  decreases  according  to  the  Indicated 
circular  radial  motion  profile  until  the  sliding  stops  at  time  t^.  We  chose  the 
circular  transition  segments  between  tg  and  t^  and  between  t^  and  t^  In  Figure 
2.2  because  the  acceleration  f  Is  proportional  to  p^  and  p2  :  we  found  this  to 
be  a  convenient  parameterization  for  studying  the  effects  of  the  radial  acce- 
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Figure  2.2:  Radial  Positioa  Profile  of  the  Sliding  >  jnlf 


One  of  the  most  Important  features  of  the  application  In  this  paper  results 
from  the  fact  that  we  specify  that  the  radial  velocity  of  the  flexible  link  stop 
very  quickly;  i.e.,  t^  -  ^2  P2  ****  small  and  the  radial  acceleration  r(t)  Is 
very  large  between  t2  and  t3.  With  this  large  i*(t).  the  term  P(t)IC3(r(t))q  m 
the  equations  of  motion  produces  the  equivalent  of  a  lateral  Impulse  on  the 
flexible  link,  unless  q(t)  Is  near  zero  between  t2  and  t3. 

Two  observations  about  (2.1)  that  are  Important  for  adaptive  control  can  be 
made  from  the  detailed  equations  of  motion  In  [Kl],  First.  (2.1)  can  be  written 
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M(t)q  +  D(t)q  +  K(t)q  ■  Bu, 


(2.7) 


where 


K(t)  -  IC(t)  -  Qg(t)B‘^g(s1n  e(t))/0(t).  (2.8) 

Second,  for  sufficiently  small  elastic  vibration  of  the  flexible  link,  no  domi- 
nant  terms  In  the  matrices  M(t),  D(t)  and  K(t)  involve  the  elastic  displacements 
q2,  q3  and  q^  or  their  time  derivatives.  Hence  the  dominant  terms  In  the  coef¬ 
ficient  matrices  In  (2.1)  and  (2.7)  vary  no  faster  than  the  rigid-body  angle 
e(t),  the  corresponding  angular  velocity,  r(t)  and  r(t).  This  conclusion 
follows  from  a  tedious  but  straightforward  examination  of  the  detailed  equations 
of  motion  In  [Kl].  It  results  essentially  from  the  fact  that  we  use  a  linearly 
elastic  model  for  the  flexible  link,  so  that  no  significant  terms  In  the 
equations  of  motion  are  more  than  first  order  In  the  elastic  degrees  of  freedom. 


Now  we  consider  digital  control  of  (2.1)  and  (2.7)  by  zero-order  sample  and 
hold;  1.e.,  at  the  beginning  of  the  k^*'  sampling  Interval  (k  •  0,  1,  2,  ...),  we 
sampl?  an  output  vector  y(k)  and  apply  a  constant  control  vector  u(k)  for  the 
duration  of  the  k^^  sampling  Interval.  We  measure  the  rigid-body  angle  and  the 
elastic  deflection  of  the  free  end  of  the  flexible  link  (the  end  holding  the 
payload),  so  that 


(2.9) 


.ccording  to  standard  linear  system  theory,  an  Input/output  model  for  (2.7) 
with  digital  Input  and  digital  linear  output  has  the  form  of  the  ARMA  model 
N  N 


y(k)  + 


A^(k)y(k-1) 


B,(k)u(k-1) 


(2.10) 
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where  A^(k)  and  B^(k)  are  matrices  of  appropriate  dimension  and  N  depends  on  the 
number  of  excited  flexible  modes.  In  our  application  A^(k)  Is  2x2  and  B^(k) 

Is  2x1.  For  our  simulations,  N  need  not  exceed  8  because  our  simulation  model 
has  four  degrees  of  freedom.  If  the  sampling  rate  is  fast  compared  to  the  time 
rates  of  change  of  the  dominant  terms  In  the  coefficient  matrices  in  (2.7)(l.e., 
If  the  sampling  rate  Is  'ast  compared  to  the  r1g1d>body  angular  velocity  and 
acceleration  and  r  and  r*),  then  the  coefficient  matrices  in  (2.10)  can  be  con¬ 
sidered  to  vary  slowly.  In  this  case,  an  adaptive  parameter  estimator  should  be 
able  to  track  the  coefficients  In  (2.10)  and  predict  y(k)  from  data  taken 
through  time  k  -  1.  Such  prediction  Is  the  basis  for  the  subsequent  adaptive 
control  algorithm. 

The  sampling  rate  used  In  the  simulations  In  Section  5  Is  100  Hz. 

Numerical  results  Indicate  that  the  slowly-varying-coefficlents  hypothesis  Is 
valid  for  our  problem  except  between  the  times  t2  and  t^  when  the  Impulsive 
effect  of  the  large  radial  acceleration  first  excites  higher  flexible  modes  that 
usually  are  not  excited,  causing  both  the  ARMA  coefficients  and  the  minimum  ARMA 
order  to  change  abruptly.  Since  t^  -  t^  Is  only  lOX  of  one  sampling  Interval, 
the  adaptive  controller  Just  sees  a  switch  from  one  set  of  slowly  varying  coef¬ 
ficients  to  another.  The  numerical  results  In  Section  5  demonstrate  that  It  is 
Important  for  the  adaptive  controller  to  be  able  to  change  both  the  ARHA  coef¬ 
ficients  and  the  ARMA  order. 

To  shed  further  light  on  the  system  to  be  controlled,  we  consider  the 
linearized  open-loop  dynamics  of  the  manipulator  with  constant  r,  payload  Mpj^  • 
.IH^,  and  no  gravity  moment  for  the  Initial  value  of  r  -  Im  and  the  final  value 
of  r  •  2m.  Table  2.2  gives  the  poles  and  zeros  of  the  continuous-time .transfer 
functions  from  the  control  torque  to  the  two  measurements.  In  the  complex  form 

i 

I 
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of  poles  and  zeros,  the  frequency  1s  given  In  rad/s;  for  the  poles,  the  damping 
ratio  C  and  the  frequency  f  In  Hertz  are  given  In  parentheses.  We  shall  empha¬ 
size  that  the  numbers  In  Table  2.2  were  not  present  explicitly  In  the  nonlinear 
finite-element  model  that  we  used  for  all  simulations.  We  computed  these  poles 
and  zeros  after  we  linearized  the  simulation  model. 

As  is  cormon  with  flexible  structures,  we  have  a  mimimum-phase  transfer 
function  for  the  sensor  colocated  with  the  actuator  (I.e.,  the  hub-rotation  sen¬ 
sor),  The  transfer  function  for  the  noncolocated  sensor,  which  measures  the 
elastic  tip  deflection^ Is  also  minimum-phase  after  the  double  pole-zero  can¬ 
cellation.  (The  transfer  function  from  the  control  torque  to  the  absolute  tip 
displacement  for  r  >  Im  has  a  real  zero  at  -t-lTS,  but  our  controller  does  not  use 
the  absolute  tip  displacement.) 

Another  common  feature  In  digital  control  of  flexible  structures  Is  that, 
whatever  the  sampMng  rate,  there  will  be  some  modes  with  frequencies  above  the 
Nyquist  frequency.  Thus  we  have  selected  the  parameters  In  Table  2.1  so  that, 
for  ■  .IM2  and  r  •  Im,  the  linearized  plant  has  two  frequencies  above  the 
Nyquist  frequency  (50  Hz  In  our  problem)  and,  for  r  ■  2m,  there  Is  one  plant 
frequency  above  the  Nyquist  frequency.  As  Illustrated  by  the  simulations 
described  In  Section  5,  the  modes  with  the  two  nighest  frequencies  have  suf¬ 
ficient  open-loop  damping  that  they  are  Important  in  the  control  problem  only 
after  being  excited  by  the  Internal  Impulse  at  t2* 

For  any  value  of  r,  the  continuous-time  linearized  system  Is  controllable 
because  there  are  no  open-loop  mode  shapes  In  which  the  hub  rotation  Is  zero  and 
there  are  no  repeated  open-loop  eigenvalues  besides  zero,  which  corresponds  to 
the  rigid-body  mode.  For  any  constant  value  of  r  at  which  no  frequency  Is  a 
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multiple  of  50  Hz,  the  discrete-time  1,:,...  a>ite.T.  (Incl the  mcces  wUh 

frequency  above  50  Hz)  Is  controllable  In  the  sense  that  the  state  vector  can  be 
driven  to  zero  1n  a  finite  number  of  steps.  With  the  positive  damping  modeled 
In  the  beam,  even  the  nonlinear  plant  Is  stablllzable  for  any  constant  value  of 
r.  These  conclusions  are  all  easy  to  verify.  We  believe  and  all  of  our  simula¬ 
tions  Indicate  that  similar  statements  about  controllability  and  stablllzablllty 
hold  for  time-varying  r,  but  we  have  not  done  the  rigorous  analysis. 

While  ve  have  eelecced  the  paraaeters  in  the  aioulation  aodel  to  produce 
plant  characteristics  that  are  conmon  In  control  problems  for  flexible 
structures,  this  paper  does  not  address  the  questions  that  arise  when  a 
sequence  of  finite  element  models  with  increasing  dimensions  is  used  for 
control  design  and  simulation.  Using  more  elements  in  the  simulation  model  in 
this  paper  would  add  more  high  frequency  modes  above  the  Nyquist  frequency, 
and  under  some  excitations  these  additional  modes  might  require  the  order  of 
the  controller  to  be  larger  than  the  maximum  order  used  in  the  simulations 
here.  The  question  of  what  order  the  finite  element  models  used  for 
controller  design  and  simulation  mvist  have  so  that  the  controller  can  be 
guaranteed  to  work  well  for  a  distributed  model  of  a  flexible  structure 
requires  approximation  theory  and  convergence  analysis  beyond  the  scope  of 
this  paper.  Such  Issues  have  been  addressed  extensively  for  optimal  control 
of  flexible  structures  (see  [GAl,  GA2]  for  example),  end  corresponding  theory 
for  adaptive  control  of  flexible  structures  is  e  focus  of  current  research. 

The  purpose  of  this  paper  is  to  demonstrate  the  desirability  of  e  controller 
that  can  vary  its  order  adaptively  to  account  for  modes  that  sometimes  ere 
excited  end  sometimes  not,  end  to  demonstrate  how  such  e  controller  can  be 
constructed. 
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TABLE  2.2 


Poles  and  Zeros  for 


Continuous-Time 
Mpj_  ■  0.1  X  M2 


r  ■  Im 

Poles 

0 

0 

-2.8  ±  1  74.7 
-137.5  ±  1506.1 
-1278.7  ±  1960.4 


Zeros 
Channel  1 


-2.3  ±  1  67.8 
-132.7  ±  1497.8 
-1263.7  ±  1964.6 


(C.O 


(.037,  11.9H2) 
(.272,  80.5H2) 
(1.33,  152. 9H2) 


Zeros 
Channel  2 

0 

0 

-173.5  ±  1563.0 
-1135.4  ±  1990.8 


r  ■  Zm 

Poles 

0 

0 

-0.4  ±  1  26.8 
-12.1  ±  1155.1 
-102.5  ±  1441.0 


Zeros 
Channel  1 


-0.3  ±  1  22.8 
-11.3  ±  1149.7 
-100.1  ±  1436.1 


(C.O 


(.015,  4.3H2) 
(.078,  24.7H2) 
(.232,  70.2H2) 


Zeros 
Channel  2 

0 

0 

-15.6  ±  1176.2 
-86.7  ±  1407.2 
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3.  PARAMETER  IDENTIFICATION 


For  adaptive  Identification  of  the  coefficient  matrices  A^(k)  and  B^(k)  In 
(2.10)  and  determination  of  the  ARMA  order  N,  we  use  the  least  squares  lattice 
algorithms  In  Tables  3.1  and  3.2.  These  algorithms  are  from  [JGl],  and  similar 
algorithms  can  be  found  In  [lMFI,  FI,  HMl]  and  other  references. 

We  will  discuss  the  structure  of  the  lattice  filter  only  briefly  to  Indi¬ 
cate  the  most  Important  points  for  the  application  here.  For  more  detail  see 
[FI,  HMl,  Jl,  JGl,  LSI].  The  lattice  structure  Is  based  on  two  sets  of  vectors 
called  forward  and  backward  residual  errors.  The  forward  residual  error  vectors 
are  obtained  from  projection  of  the  most  recent  regression  vectors  onto  the  span 
of  previous  regression  vectors.  (Regression  vectors  contain  measurement 
histories.  In  the  lattice  filter,  both  control  system  Inputs  and  outputs  must 
be  treated  as  measurements.)  The  order  of  the  lattice  Is  the  order  of  the  ARMA 
model  with  which  the  lattice  attempts  to  fit  the  data.  The  norm  of  the 
order  forward  residual  error  Is  the  minimum  value  of  the  objective  functional 
to  be  minimized  by  a  least-squares  estimate  of  the  parameters  In  an  ARMA  model 
of  order  N.  The  backward  errors  are  a  set  of  Gram-Schmidt  vectors  that  span  the 
same  space  as  the  regression  vectors.  At  each  step,  the  filter  updates  the 
first  element  of  each  of  the  error  vectors  and  the  norms  of  the  error  vectors. 
The  N^^-order  lattice  generates  the  equivalent  of  the  least-squares  fit  to  data 
for  every  ARMA  order  between  1  and  N. 

The  main  lattice  algorithm,  to  be  run  at  every  time  step.  Is  listed  In 
Table  3.1.  The  1x3  matrices  d|^(k)  and  r|^(lc)  are  the  first  components  of  the 
forward  and  backward  residual  errors,  respectively.  See  [JGl].  Also,  Rjj(k), 
R{j(k),  K^(k)  are  3x3  matrices  and  G|^(k)  Is  a  scalar.  The  forgetting  factor  A  Is 
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a  positive  real  number  less  than  or  equal  to  1,  which  (when  A  <  1)  reduces  expo¬ 
nentially  the  effect  of  older  data.  For  the  simulations  In  Section  5,  we  used 

A  ■  *99. 


The  statement  that  the  lattice  Is  order  recursive  refers  to  the  fact  that 
the  maximum  lattice  (ARMA)  order  can  be  Increased  by  1  at  each  time  step,  up  to 
some  limit  determined  by  on-line  computing  capacity.  In  practice.  N  is 
Increased  until  either  It  reaches  the  upper  limit  or  some  criterion  Involving 
the  matrix  R^(k)  Indicates  that  N  need  not  be  Increased  further. 

The  diagonal  elements  of  the  matrix  R^Ck)  are  the  squares  of  the  norms  of 
the  forward  residual  errors,  which  are  related  closely  to  prediction  error  (see 
[JGl]).  These  diagonal  elements  Indicate  the  degree  to  which  an  N^^-order  ARHA 
model  fits  the  data,  and  they  can  be  used  to  determine  the  order  of  the  plant. 
In  particular,  the  (1,1)  element  of  R^(k)  Is  a  measure  of  the  accuracy  of  the 
channel  of  the  ARMA  mode!  corresponding  to  the  rigid-body  angle  of  our  manipula 
tor.  We  will  use  this  number  to  determine  the  order  of  the  ARMA  model  to  be 
used  for  adaptive  control,  as  discussed  in  Section  4.5. 


The  ARMA  coefficients  are  obtained  from  the  algorithm  In  Table  3.2,  where 
\  l(l«)f  ^(k)  and  ^(k)  are  3x3  matrices.  For  each  N,  the  matrix  Ajj  ^(k) 
has  the  form 


-  A)(fc)^l  * 

B,(k)^|  X 


(3.1) 


Where  A^(k)  and  B^(k)  are  the  least-squares  estimates  at  time  k  of  the  matrices 
In  (2.10). 
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Initialize; 


/Jo'(-l)=^-n  =  <T/.  Ks^m=0  for  A'+l>i^O 
where  <t  is  a  small  number  (e.g.,  a  =  10’^) 

For  each  A  ^  0: 

e^k)  =  ro(A:)  =  ly\k)  ii(A:)] ,  C(f,k)  =  I 
/?o'W=  P^k)=edikf  ef/ik)+}.P^k-  1) 

For  each  k^  I,  for  0  to  k-  I: 

^v+i(*)=  J)+  iy,,(A:-  1) 

ef,k)-r^k-  -  I)A'^+,(A:) 

''/v+i(^)  =  ''/s<*-  1)- 
/?v+iW= 

RU i<*)  =  -  I)  -  ,(^)^.v'(^V^/v+ 1(*) 

C/v+i(*)=  r;,(A)/?;/(A:y^(A:) 

Table  3.1  Residual  Error  Lattice  Algorithm 


For  1,  A:  -  1  and  for  f=  I,  N: 

Cn*  I  ;(*) »  Cyv/*)  -  B^fk)RJ,\kyll,k) 

1 , f+ i(Ac)  =  C’/v/^)^^*(A;y^A:)]  —  Afffk)Rff  (Ac)/f^+j(Ac) 

with 

RN*^k)Kf/^f.k) 

cWi.^v+iW*  Ri/ikyjjf^k) 

Table  3.2  Algorithm  For  AR  Coeiricienis 


The  adaptive  controller  in  this  paper  is  cooputationally  efficient 
because  the  bulk  of  the  computation  is  in  the  lattice  filter  and  lattice 
filters  are  the  fastest  digital  signal  processing  algorithms  for  recvirsive 
least-squares  estimation.  The  signal-flow  structure  of  the  lattice  lends 
itself  particularly  well  to  parallel  architectures.  These  points  are 
discussed  widely  in  the  digital  signal  processing  literature  (e.g.,  [FI,  GSl, 
HMl,  LMFl,  LSI]). 
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4.  ADAPTIVE  CONTROL 


4.1  Control  Law 


The  adaptive  control  law  Is  based  on  (2.10)  with  estimated  coefficient 

A  A 

matrices  A^(k)  and  B^(k).  For  defining  and  computing  the  control  law,  the  / 


ARMA 


order  N  is  assumed  to  be  fixed;  changing  N  adaptively  Is  discussed  In  Section 


We  set 

N  N 

c(k)  -  Y-  A.(k)  y(k-1)  -  H  B.(k)  u{k-1)  -  y  y{k-l)  +  (1+Y)y  (k)  (4.: 

1-1  ’  1-1  ^ 

where  Y  Is  a  real  number  with  r..jgn1tude  less  than  1  and  y^(k)  Is  a  reference 
signal  to  be  discussed  later.  We  choose  u(k>l)  to  minimize 


J(k)  .  e(k)V(k)  +  nou2(k-l) 


Where  Q  Is  the  nonnegative  matrix 


(4.2) 


(4.3) 


and  Hq  Is  a  nonegative  real  number.  When  (Hq  B^(k)  QB^(k))  >  0,  the  unique 
control  u(k-l)  that  minimizes  (4.Z)  Is 

u(k-l)  . 

N  N 

B}(k)0{Z!  >(k-1)  -  H  B^(k)  u{t-1)  -  Y  y(k-l)  ♦ 


Hg  +  Bj(k)0Bj(k) 


(4.4) 
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The  ODjectIve  J(k)  Is  motivated  by  the  fact  that  1f  €(k)  ■  0,  then  (4.1) 
and  (2.10)  yield 

y(k)  +  Y  y(k-l)  •  (1+Y)yp(k).  (4.5) 

Since  In  our  problem  y(k)  and  yp(k)  are  2-vector5  and  u(k-l)  Is  a  scalar,  it  is 
almost  never  possible  to  obtain  e(k}  ■  0.  On  the  other  hand,  If  >  0  and  ri2  * 
Hq  ■  0,  the  first  element  of  c(k)  will  be  zero,  so  that  (4.5)  will  hold  with 
y(k)  and  yj.(k)  replaced  by  their  respective  first  elements.  (The  first  element 
of  is  nonzero  in  our  problem.) 

The  adaptive  control  law  In  (4.4)  Is  a  variation  on  a  family  of  model 
reference  schemes  discussed  In  [GSl,  Section  6.3].  For  a  plant  that  can  be 
represented  exactly  by  (2.10)  with  fixed  order  N  and  constant  coefficients,  sta¬ 
bility  results  for  the  closed-loop  system  produced  by  the  adaptive  controller 
here  are  similar  tc  stability  results  in  [GSl,  Chapters  5  and  6].  In  par¬ 
ticular,  1f  the  first  element  of  Bj  Is  nonzero,  tIj  >  0  and  112  ■  Hq  ■  0,  then  the 
adaptive  controller  here  reduces  to  a  one-step-ahead  model  reference  adaptive 
controller  for  a  SISO  system.  In  that  case,  a  sufficient  condition  for  asymp¬ 
totic  stability  Is  that  all  plant  zeros  lie  Inside  the  open  unit  disc. 

While  (2.10)  with  constant  coefficients  cannot  represent  the  manipulator  In 
our  problem  exactly,  we  have  considered  linearized  motion  about  the  final 
equilibrium  position  for  asymptotic  stability  analysis.  (In  this  case,  the 
flexible  link  no  longer  slides,  and  the  control  torque  Is  perturbed  about  the 
appropriate  static  torque.)  A  root  locus  analysis  in  [Kl]  shows  how  the 
(dlscictf'tlme)  zeros  of  the  open-loop  ARMA  model  that  relates  torque  pertur¬ 
bations  to  perturbations  in  the  rigid-body  angle  depend  on  the  viscoelastic 
damping  In  the  flexible  link.  All  of  these  zeros  lie  on  the  unit  circle  when  no 
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damping  is  modeled  In  the  open-loop  manipulator,  and  all  of  these  2eros  move 
inside  the  unit  circle  when  damping  Is  modeled  In  the  flexible  link.  This  root 
locus  analysis  Is  straight  forward.  The  analogous  distributions  of  continuous¬ 
time  zeros  for  flexible  structures  with  colocated  actuators  and  sensors  is  well 
known. 


For  controlling  the  manipulator,  we  take  ■  1000  and  112  ”  1*  Hence  the 
control  law  In  (4.4)  almost  amounts  to  SISO  control  of  the  rigid-body  angle, 
except  near  the  final  position,  where  the  relatively  small  penalty  on  the 

_Q 

elastic  tip  deflection  slightly  Improves  settling.  Also,  we  take  Hq  ■  SxlO  in 
(4.2)  and  (4.4)  so  that  the  control  law  will  be  determined  uniquely,  and  because 

_Q 

a  root  locus  analysis  In  [Kl]  Indicates  that  values  of  tIq  on  the  order  of  10  , 
as  opposed  to  Hq  "  0,  Improve  the  asymptotic  stability  of  the  closed-loop  system 
without  producing  measurable  steady-state  error. 

During  the  first  10  steps  when  the  manipulator  starts  to  rotate  (I.e., 
after  the  learning  period),  we  set  Y  •  0  to  get  the  manipulator  moving  rapidly. 
Then  we  set  v  ■  0.96  for  the  remainder  of  the  motion  so  that  the  controller  does 
not  try  to  eliminate  the  error  between  y(k)  and  y^(k)  unrealistically  fast. 

4,2  Reference  Signal 

We  are  most  Interested  In  commanding  the  absolute  position  of  the  manipula¬ 
tor  tip,  which  lies  on  a  circle  of  radius  L  -f  r.  Since  r(k)(k  •  0,1,2,  ...)  is 
already  specified  by  Figure  2.2,  we  concern  ourselves  with  the  tip  position 
ytip(k)  measured  as  arc  length  along  the  circle  of  radius  L  4  r(k).  Since 

nipt'')  -  e(k)(L4r(k))  4y2(k).  (4.6) 

we  command  y^^p(k) (whose  desired  profile  Is  shown  in  Figures  5. 1-5.6)  indirectly 
through  the  reference  signal  yp(k)  for  the  output  vector  in  (2.9). 
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In  8p(t),  the  commanded  rigid-body  angle,  we  must  account  for  the  static 
tip  deflection  due  to  gravity.  We  cannot  determine  this  deflection  off-line 
without  using  the  exact  data  for  the  flexible  link  and  the  exact  payload. 
Therefore,  to  give  the  control  system  the  capability  to  adapt  to  unknown  link 
characteristics,  different  payloads  and  different  desired  final  positions,  we 
use  an  estimate  of  the  static  tip  deflection  in  the  reference  signal.  We 
construct  this  estimate,  denoted  by  low-pass  filter 

y2(k+l)  +  ay2(k)  •  (l+a)v(k), 

v(k)  -  Slgn(y2(k))*m1n 

which  attenuates  oscillations  In  y2(l()*  The  constants  In  (4.7)  are  a  ■  -0.986 
and  ■  0.1  m.  During  the  large-angle  motion.  It  1s  Important  only  that  y2(l() 
be  small;  near  the  final  position  should  approximate  closely  the  steady- 

state  tip  deflection. 

We  denote  the  polar  coordinates  of  the  desired  final  tip  position  by  (0^, 
L^)  where  ■  L  +  r^  and  r^  Is  the  final  value  of  r.  The  reference  signal  for 
the  output  vector  In  (2.9)  Is 


with 

e^Ck)  -  -  y2(k)]/(Ur(k)).  (4.9) 

This  definition  of  6^(k)  Is  motivated  by  and  equivalent  to 

.  9f(k)(Ur(k))  .ijCk),  .  (4.10) 

Which  Is  obtained  from  (4.6)  by  replacing  y2(k)  with  the  estimate  y2(k)  and 
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replacing  desired  tip  position 

We  acknowledge  that  the  reference  signal  here  is  unusual  because  It  in¬ 
cludes  Implicitly  a  sensor  measurement.  We  see  no  way  to  avoid  this  feature  if 
the  control  system  Is  to  be  truly  adaptive. 

4.3  Learning  Period 

To  allow  the  lattice  filter  to  obtain  Initial  parameter  estimates  before 
the  manipulator  starts  to  move  significantly,  we  use  an  Initial  learning  period 
of  0.2  sec.  or  20  samples.  During  this  period,  the  control  torque  Is 

Uiearn^*')  *  cos(.4ffk)  Nm.  (4.11) 

Because  of  the  frequency  of  this  Input,  the  manipulator  does  not  move  signifi¬ 
cantly  from  Its  Initial  position  at  6  ■  0.  We  do  not  attempt  to  Identify  the 
plant  order  or  parameters  precisely  during  this  learning  period  because  the  ARHA 
model  will  begin  to  change  as  the  flexible  link  begins  to  slide  and  nonlinear 
terms  Involving  the  angular  velocity  build  up  during  the  fast  rigid-body  rota¬ 
tion.  Since  the  flexible  modes  are  not  excited  significantly  when  the  manipula¬ 
tor  starts  to  rotate  at  .2  sec,  we  use  the  ARHA  order  2  during  and  Immediately 
after  the  learning  period.  Therefore,  we  need  only  a  simple  Input  for  a  short 
learning  period  to  obtain  parameter  estimates  sufficient  for  beginning  the 
rigid-body  motion. 

4.4  Bound  on  Control  Magnitude 

In  the  simulations  In  Section  5,  we  Impose  the  bound  u^^^  •  5x10^  Mn  on  the 
magnitude  of  the  control  torque  to  demonstrate  that  the  algorithm  Is  robust  with 
respect  to  actuator  saturation.  This  means  that  if  the  magnitude  of  the  control 
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torque  In  (4.4)  exceeds  then  times  the  sign  of  the  expression  on  the 

right  side  of  (4.4)  is  used  for  u(k-l). 

4 

The  torque  limit  5x10  Nm  Is  approximately  the  torque  required  to  produce 
the  angular  acceleration  associated  with  the  reference  signal  In  the  simula¬ 
tions.  It  Is  not  clear  whether  such  a  torque  limit  is  realistic  for  robot 
actuators;  however,  our  objective  In  the  simulations  is  to  demonstrate  an  adap¬ 
tive  control  algorithm  that  can  control  fast  slewing  in  the  presence  of  elastic 
vibrations  and  a  variable-order  plant.  If  lower  torque  limits  on  available 
actuators  dictate  slower  angular  accelerations,  the  adaptive  control  problem 
should  be  less  challenging  because  the  slower  rigid-body  motion  will  excite  the 
elastic  modes  less  than  In  our  simulations. 

4.5  Changing  the  Order  of  the  ARMA  Model 

As  we  said  In  Section  3,  the  size  of  the  diagonal  elements  of  R^(k} (which 
must  be  nonnegative)  Indicates  how  well  an  ARMA  model  of  order  N  fits  the  Input- 
output  data  taken  through  time  k.  In  [JGl],  data  from  a  simulated  flexible 
structure  wa$  used  to  demonstate  how  R^(k)  Indicates  the  number  of  excited  modes 
when  small  measurement  noise  Is  present.  The  idea  is  to  increase  N  and  look  for 
large  drops  In  the  diagonal  elements  of  R^(k).  It  might  appear  then  that  R^(k) 
should  be  examined  to  determine  the  plant  order  detected  by  the  lattice  filter, 
and  that  that  ARMA  order  should  be  used  for  adaptive  control.  However,  order 
determination  for  adaptive  control  Is  not  so  simple,  for  two  reasons.  First,  an 
ARMA  model  with  constant  coefficients  cannot  fit  exactly  the  nonlinear,  time- 
varying  manipulator  to  be  controlled  here  (even  without  neasurement  noise),  so 
that  we  should  not  expect  the  kind  of  sharp  drop  In  R^Ck)  at  the  correct  order 
that  was  seen  In  [JGl].  Second,  as  demonstrated  In  [JG2]  with  data  from  an 
experimental  flexible  structure,  the  lattice  filter  can  detect  very  marginally 
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excited  flexible  modes,  which  usually  do  not  need  to  be  controlled  actively. 

Our  simulations  have  Indicated  repeatedly  that  using  an  ARHA  order  that 
corresponds  to  the  total  number  of  modes  in  the  plant  leads  to  poor  iden¬ 
tification  and  control  1f  some  of  the  modes  are  only  marginally  excited.  While 
marginally  excited  modes  were  identified  In  [JGl,  JG2],  much  longer  data  records 
were  required  than  an  adaptive  controller  has  time  to  process  before  the  control 
law  must  be  computed  and  executed. 

Our  criterion  for  choosing  the  ARMA  (lattice)  order  adaptively  after  the 
learning  period  combines  information  In  R^(k)  with  a  measure  of  the  control 
system  performance  to  obtain  an  Indication  of  the  order  of  the  ARHA  model  needed 
for  effective  control.  The  (1,1)  element  of  R^(lc),  denoted  by  R2(k,l,l),  Is  the 
term  In  R^(k)  most  closely  related  to  the  accuracy  of  one-step-ahead  prediction 
for  the  rigid-body  measurement  In  our  problem  (again,  see  [JGl]),  and  most  per¬ 
tinent  for  our  controller  because  of  the  relative  weighting  In  the  matrix  Q.  At 
each  sampling  time,  the  following  test  Indicates  whether  to  change  N: 

e(k)  -  - £ -  ;  (4.12) 

k 

IAq  +  ^  |y(J)  -  y.(J)l^J 

^  j-0  ^ 

R5(k,l,l)  <  e(k)  -  A.  — >  N  -  N  -  1. 

"  *  (4.13) 

Rj(k,l.l)  >  e(k)  ♦  Aj  — >  N  -  N  4  1. 

-S 

We  have  chosen  the  constants  c  •  50.0,  Aq  •  1.0x10  ,  and  A^  •  0.005 
empirically. 

This  order-change  criterion  Is  admittedly  ad  hoc,  and  further  research 
might  produce  a  more  sophisticated  criterion.  However,  the  test  here  has  two 
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features  that  we  believe  will  be  Important  In  any  criterion  for  changing  the 
order  of  an  adaptive  controller;  1)  the  model-fit-to-data  error  can  remain  large 
If  the  control  system  achieves  the  desired  response;  2)  the  dead  band  repre¬ 
sented  by  Aj  reduces  chattering  of  the  ARMA  order  significantly  (although  the 
simulations  In  Section  5  show  that  order  chattering  has  not  been  eliminated 
entirely).  Because  magnitudes  of  desired  and  true  responses  and  number  of 
sampling  points  for  a  typical  motion  vary  with  the  application,  we  suspect  that 
any  criterion  with  features  1)  and  2)  will  have  constants  that  must  be  adjusted 
empirically  for  each  application. 

5.  SIMULATION 

In  the  simulations,  the  nonlinear  equations  of  motion  (2.1)  were  solved 
numerically  with  the  control  torque  generated  by  the  adaptive  control  law  In 
(4.4)  and  the  reference  signal  given  In  Section  4.2.  (See  [Kl]  for  details  of 
the  numerical  Integration).  The  parameters  In  the  ARMA  model  (2.10)  were  esti¬ 
mated  with  the  lattice  algorithms  In  Tables  3.1  and  3.2,  and  the  order  of  the 
adaptive  controller  (I.e.,  the  order  of  the  ARMA  model  upon  which  the  control 
law  Is  based)  was  determined  adaptively  by  the  criterion  In  Section  4.5.  Also, 
the  learning  period  and  control  torque  bounds  1n  Sections  4.3  and  4.4  were  used. 

Representative  examples  of  the  numerous  simulations  In  [Kl]  are  presented 
In  Figures  5. 1-5.6.  For  these  figures,  the  initial  position  Is  the  vertical 
position  corresponding  to  6  ■  0  In  Figure  2.1.  In  the  final  position. 
Illustrated  by  the  all-dotted  position  In  Figure  2.1,  the  polar  coordinates  of 
the  manipulator  tip  are  6^  •  2  rad  and  •  5m  (recall  Section  4.2).  In  the 
figures.  TIP  POSITION  is  the  y^^p  In  (4.6)  and  DESIRED  TIP  POSITION  at  time  k  is 
e^L^(l-e'*°®*‘) (recall  (4.9)  and  (4.10)).  TIP  POSITION  Is  In  meters,  and  TIME  is 
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1n  seconds.  CONTROL  TORQUE  In  the  figures  Is  the  control  u  divided  by  the  bound 


For  all  six  figures  here,  the  same  adaptive  controller  was  used  (i.e.,  all 
the  same  constants  In  the  equations  In  Section  4),  except  that  In  Figures  5.2 
and  5.4,  the  order  N  of  the  adaptive  controller  was  fixed  at  2.  For  Figures 
5.1,  5.3,  5.5,  and  5.6,  the  maximum  allowable  order  of  N  was  8.  To  emphasize 
the  role  that  the  lattice  filter  plays  In  adjusting  the  order  of  the  controller, 
the  order  of  the  controller  Is  labeled  LATTICE  ORDER  In  the  figures. 

The  gravitational  moment  on  the  manipulator  was  Included  In  the  simulation 
model  for  Figures  5. 1-5.4,  but  no  gravity  was  modeled  for  Figures  5.5  and  5.6. 

In  other  words,  the  manipulator  moved  In  a  vertical  plane  for  Figures  5. 1-5.4 
and  In  a  horizontal  plane  for  Figures  5.5  and  5.6.  Thus  a  steady-state  control 
torque  equal  to  approximately  1%  of  the  maximum  control  torque  Is  required  In 
Figures  5. 1-5.4  to  offset  the  torque  due  to  gravity.  The  steady-state  control 
torque  In  Figures  5.5  and  5.6  Is  zero. 

Figures  5.1  and  5.2  compare  the  variable-order  adaptive  controller  with  an 
adaptive  controller  of  fixed  order  2  for  a  payload  equal  to  lOX  of  the  mass 
of  the  flexible  link.  Until  the  effective  lateral  Impulse  at  0.8  sec  produced 
by  the  sudden  halt  of  the  radial  motion,  the  second-order  adaptive  controller  is 
sufficient.  The  Impulse  greatly  excites  transient  elastic  vibration  In  the 
manipulator  and  the  lattice  (ARHA)  order  In  Figure  5.1  Is  adjusted  automatically 
(In  six  steps)  to  the  maximum  order  N  «  8.  When  most  of  the  vibration  Is  taken 
out  of  the  system,  the  lattice  order  Is  decreased  back  to  2. 

It  Is  not  surprising  that  a  second-order  controller  produces  poor  transient 
response  In  the  presence  of  significant  vibrations  of  the  flexible  link.  A  more 
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Interesting  question  Is  whether  we  should  not  Just  fix  the  order  of  the 
controller  at  8  throughout  the  motion.  We  have  simulated  the  response  for 
fixed-order  controllers  with  N  ■  8  and  other  orders  between  2  and  8,  and  the 
simulations  almost  always  show  a  response  that  becomes  unstable  before  .8  sec, 
when  the  higher  order  Is  needed  In  the  controller.  It  appears  that  when  the 
order  of  the  ARMA  model  Is  higher  than  that  needed  to  fit  the  Input/output  data 
for  the  plant,  the  estimates  of  the  redundant  parameters  are  so  poor  that  the 
adaptive  prediction  upon  which  the  control  law  Is  based  Is  very  poor. 

Figures  5.3  and  5.4  compare  the  variable  and  fixed-order  adaptive 
controllers  for  a  payload  equal  to  50%  of  the  mass  of  the  flexible  link.  While 
the  superiority  of  the  variable-order  controller  is  still  clear  from  Figures  5.3 
and  5.4,  these  and  other  simulations  in  [Kl]  show  that  the  difference  between 
the  performance  produced  by  the  variable-order  controller  and  that  produced  by 
the  fixed-order  controller  decreases  with  Increasing  payload.  With  a  payload 
equal  to  70%  of  the  mass  of  the  flexible  link,  both  the  variable  and  the  fixed- 
order  controllers  produce  a  response  very  similar  to  that  In  Figure  5.3  (when 
the  gravitational  moment  Is  present). 

We  are  not  certain  why  the  change  In  controller  order  1$  less  Important  for 
heavier  payloads.  The  need  to  Increase  the  controller  order  Is  greatest  when 
the  Internal  Impulse  at  0.8  sec  has  Its  greatest  effect  on  the  transverse  vibra¬ 
tions  of  the  flexible  link,  and  the  magnitude  of  this  Impulse  Is  proportional 
to  the  magnitude  of  the  elastic  deformation  at  the  time  of  the  Impulse.  More 
extensive  data  In  [Kl]  Indicates  that,  for  the  commanded  motion  In  the  simula¬ 
tions  here,  the  elastic  deformation  with  the  larger  payloads  Is  small  at  0.8 
sec.  For  heavier  payloads,  when  we  time  the  Internal  Impulse  to  coincide  with 
larger  elastic  deformation,  the  Improvement  made  by  the  variable-order  control- 
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ler  over  the  fixed-order  controller  becomes  larger;  however.  In  most  simulations 
the  difference  between  variable-order  and  fixed-order  controllers  Is  still 
greater  for  lighter  payloads. 

The  control  constants  Hqi  <12*  ^tc.  In  the  adaptive  control  law  were 
selected  empirically  to  produce  good  response  for  a  wide  range  of  payloads  when 
the  manipulator  moves  In  a  vertical  plane  under  the  Influence  of  gravity. 

Figures  5.5  and  5.6  demonstrate  that  the  variable-order  adaptive  controller  Is 
sufficiently  robust  to  continue  to  stabilize  the  manipulator  about  the  desired 
final  position  when  the  gravitational  torque  Is  removed  from  the  simulation 
model,  although  the  response  does  degrade.  The  control  constants  can  be 
adjusted  to  optimize  the  response  In  the  horizontal  plane,  but  using  the  same 
controller  for  all  of  the  Figures  here  better  demonstrates  the  adaptive  capabi¬ 
lity  of  the  controller. 

Figures  5.1  and  5.5  show  some  chatter  In  the  lattice  order.  The  criterion 
In  Section  4.5  for  order-determination  has  eliminated  most  such  chatter.  A 
better  criterion  might  eliminate  the  chatter  entirely.  We  have  no  simulations 
with  the  order  determination  criterion  here  In  which  this  chatter  appears  to 
degrade  the  response  of  the  manipulator. 

6.  CONCLUSIONS 

The  motion  of  the  manipulator.  In  both  the  horizontal  and  vertical  planes, 
can  be  controlled  adaptively  for  a  wide  range  of  payloads.  Because  of  the 
Internal  Impulse  associated  with  the  sudden  halt  In  sliding  of  the  flexible 
link,  the  capability  of  the  controller  to  vary  its  order  adaptively  Is  very 
important.  The  lattice  filter  used  for  the  on-line  parameter  estimation  pro- 
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vides  the  variable-order  capability,  which  the  standard  least-squares  algorithm 
and  other  parameter  estimation  schemes  commonly  used  In  adaptive  control  do  not 
provide.  Also,  since  a  lattice  of  order  N  generates  parameter  estimates  for 
ARMA  models  of  all  orders  between  1  and  N,  the  additional  computation  required 
for  the  controller  capable  of  varying  its  order  between  1  and  N,  as  opposed  to  a 
controller  of  fixed  order  N,  Is  only  the  minor  computation  required  for  the 
order  determination  criterion  In  Section  4.5. 
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ABSTRACT 

For  linear  systems  involving  uncertain  parameters  with  known, 
constant  nominal  values  and  uncertain  perturbations  that  vary 
sinusoidally  with  time,  Lyaptmov  robustness  analysis  is  used  to  determine 
a  stability  bound,  or  margin,  for  the  amplitudes  of  the  parameter 
perturbations.  This  bound  is  the  size  of  a  hypercube  in  parameter  space 
for  which  asymptotic  stability  is  guaranteed.  The  bound,  which  is  based 
on  a  quadratic  Lyapunov  function  that  depends  linearly  on  parameter 
perturbations,  varies  with  the  frequency  of  the  uncertain  parameter 
perturbations.  The  bound  is  asymptotically  proportional  to  the  square 
root  of  this  frequency  as  it  becomes  large. 
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1.  Introduction 

Numerous  recent  papers  [Bl-KBHl,  PT1*ZK2]  have  used  quadratic  Lyapunov 
functions  to  develop  robustness  bounds  for  linear  systems  with  uncertain 
parameters.  Some  papers  [Hfi2,  PTl,  Yl-ZKl]  have  dealt  with  robustness 
analysis  only,  while  some  [Bl.  HBl,  KBHl,  FI,  PHI,  ZK2]  have  used  Lyapunov- 
based  robustness  analysis  as  a  basis  for  design  of  robust  controllers.  A 
common  feature  of  the  references  Just  cited  and  most  related  work  is  chat  a 
single  Lyapunov  function  is  tised  for  the  entire  set  of  parameters  for  which 
stability  is  guaranteed.  Because  of  this,  the  Lyapunov  robustness  analysis 
applies  Co  time-varying  uncertain  parameters  (although  Che  nominal  plant  must 
be  constant).  However,  Che  robustness  bounds,  or  margins,  produced  by  such 
analysis  involve  only  the  magnitude  of  parameter  variations;  the  analysis 
cannot  detect  how  the  allowable  magnitude  of  \incercain  time -varying  parameters 
depends  on  their  frequency. 

In  [LI,  LGl] ,  a  quadratic  Lyapunov  function  was  developed  that  varies 
linearly  with  uncertain  plant  parameters.  Because  of  the  linear  dependence  on 
parameters,  the  method  in  [LI,  LGl]  is  called  a  first-order  method.  For  all 
but  one  example  to  date,  this  first-order  method  has  yielded  larger  robustness 
bounds  Chan  the  sharpest  possible  method  based  on  parameter- independent 
Lyapunov  functions  (see  [LI.  LGl]).  The  first-order  method  in  [LI,  LGl]  does 
not  apply  Co  problems  with  time-varying  uncertainties,  though- 

This  paper  extends  the  approach  in  (LI,  LGl]  to  linear  systems  in  which 
the  nominal  system  is  time -invariant  but  the  perturbations  in  uncertain 
coefficients  vary  sinusoidally  with  cine.  As  in  [U.,  LGl],  the  Lyaptmov 
function  here  varies  linearly  with  uncertain  parameters.  The  first-order  term 
in  the  Lyaptmov  matrix  satisfies  a  differential  equation  in  which  the  forcing 
tern  contains  the  sinusoidal  perturbations  from  the  nominal  plant.  As  a 
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result,  the  Lyapunov  function  and  the  resulting  robustness  margin  depend  on 
the  frequency  of  the  parameter  perturbations.  In  general,  the  robustness 
margin  is  proportional  to  the  square  root  of  the  frequency  of  the 
perturbations  at  large  frequencies. 


2.  Preliminaries 

Ve  consider  the  system 

(2.1)  i(t)  -  A(t)  x(t) 

where  x(t)  is  a  real  n-vector  and  the  n  x  n  matrix  A(t)  has  the  form 

(2.2)  A(t)  -  A(t,p)  -  Aq  +  G(p)  sinwt 

where  the  real  matrix  G(p)  is  a  linear  function  of  the  constant  parameter 
vector  p  -  (p2^  P2  ...  Pj,]^  <  and  the  real  matrix  Aq  is  constant  and 
independent  of  p. 

Definition  2,1.  A  real  n  x  n  matrix  function  P(t)  is  a  Lyspunov  matrix  for 
A(t)  if  i)  F(t)  is  periodic  with  period  2n/u;  il)  P(t)  is  symmetric  and 
positive  definite  for  each  t,  its  maximum  eigenvalue  is  bounded  uniformly  in  t 
and  its  minimum  eigenvalue  is  bounded  away  from  0  uniformly  in  t;  iii)  P(t)  is 
piecewise  continuously  differentiable  and  the  real  symmetric  matrix 

(2.3)  Q(t)  -  .(P(c)  +  A(t)’^  P(t)  +  P(t)  A(t)) 

is  nonnegative.  Furthermore,  P(t)  is  a  strict  Lyapunov  matrix  for  A(t)  if 
Q(t)  is  positive  definite  with  its  minimum  eigenvalue  bounded  away  from  0 
uniformly  in  t. 


Theorem  2.2.  (Standard  Result)  The  system  (2.1)  is  uniformly  exponentially 
stable  if  and  only  if  there  exists  a  strict  Lyapunov  matrix  for  A(t). 
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Ve  assume  that  the  eigenvalues  of  Aq  all  have  negative  real  parts,  so 
that  for  each  positive  definite  symmetric  real  n  x  n  matrix  Qq  there  exists  a 
xinique  positive  definite  symmetric  real  n  x  n  matrix  Pq  satisfying 

(2.4)  Aq  Pq  +  Pq  Aq  -  -Qq. 

Ve  will  factor  Qq  uniquely  as 

(2.5)  Qq  -  LL^ 

where  L  is  a  real  n  x  n  lower  triangular  matrix  with  positive  diagonal 
elements . 


3.  The  First-order  Method 

We  define 

(3.1)  P(t,p)  -  Pq  +  Pi(t.p) 

where,  for  each  value  of  the  parameter  vector  p,  Pj^(t,p)  is  the  unique 
periodic  solution  to 

(3.2)  ^i(t,p)  +  aJ  Pi(t,p)  +  Pi(C,p)  Aq  -  -(G^(P)  Pq  *  ^0 

That  (3.2)  has  exactly  one  periodic  solution  follows  from  the  fact  that  the 
eigenvalues  of  Aq  all  have  negative  real  parts.  The  matrix  P^Ct.p)  is  a 
linear  function  of  the  parameter  vector  p,  since  G(p)  is.  Furthermore, 
P^(c,p)  has  the  form 


(3.3) 


Pj^(t,p)  -  P^(p)  coswt  +  P^CP)  »in“t 


u 


where  the  real  symmetric  nxn  matrices  P^Cp)  and  P^j(p)  are  the  unique  nxn 
matrices  that  solve  the  equations 

(3.4) ,  -«P^(p)  +  aJp^Cp)  +  Pb(p)Ao  -  -(g'’^(p)  Pq  +  Pq  G(p)] 

(3.5)  «Pb(p)  +  aJp^(p)  +  Pa(p)Ao  -  0. 

To  motivate  our  terminology,  we  note  that,  if  the  solution  P(t.p)  to  (2.3)  for 
fixed  Q  -•  Qq  is  expanded  as  a  Taylor  series  in  p,  the  zero-order  term  is  Pq 
and  the  first-order  term  is  P2^(t,p). 

Next,  ve  define  some  quantities  that  will  be  useful  in  determining 
whether  P(t,p)  is  a  Lyapunov  matrix  for  A(t,p).  First, 

(3.6)  W(t,p)  -  L-^(G^(p)  Pi(t,p)  +  Pi(t.p)  G(p))  L-'^  siiwt 
where  L  is  the  matrix  in  (2.5).  For  a  matrix  M; 

(3.7)  singular  value  of  M; 

(3.8)  (Ml-  is  an  eigenvalue  of  M)  (for  square  M) . 


Since  Pq  Is  a  strict  Lyapunov  matrix  for  Aq,  the  following  two 
conditions,  together,  are  sufficient  for  P(t,p)  to  be  a  strict  Lyapunov  matrix 
for  A(t,p)  for  a  given  p: 

Condition  3.1.  A^y.»(Pn~^Pi  (t.p) )  is  bounded  strictly  below  1  uniformly  in  t. 
Condition  3.2.  g„.,(W(t,p))  is  bounded  strictly  below  1  uniformly  in  t. 

Definition  3 . 3 .  (Hvpercube  in  R“)  .  For  s  >  0, 

C(«)  -  (P  -  (Pi  P2  •••  PbI^*  IPil  ^  *)• 

Ue  note  that  C(s)  is  the  convex  hull  of  its  2^  vertices. 
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Definition  3.4.  If  f  Is  a  real-valued  function  defined  on  C(l),  then 
Pj^(f)  -  max  (f(p):  p  is  a  vertex  of  C(l)). 

Lemna  3.5.  Let  ^2'  ®  finite  collection  of  points  in  a  linear 

space*,  let  S  be  the  convex  hull  of  ^2»  •••  ^k^  *  f  be  a  convex 

function  defined  on  S.  Then  max  ^eS)  -  f(|j)  for  some  j. 

The  proof  of  this  lemma  is  elementary.  See  (LI]. 

Ue  recall  that  g^^yC*)  Is  a  norm  for  any  space  of  finite  dimensional 
matrices.  Hence  is  a  convex  function  on  any  such  space.  Also,  for  a 

fixed  matrix  M,  *  convex  function  on  the  space  of 

symmetric  matrices  N  of  a  given  dimension.  We  will  use  these  facts,  along 
with  Lemma  3.5,  to  estimate  the  largest  hypercube  C(s)  such  that,  for  each  p 
in  the  interior  of  C(s),  P(t,p)  is  a  strict  Lyapunov  matrix  for  A(t,p). 

Since  P^(p)  and  P^Cp)  sre  linear  in  p  and  since  a  convex  function  of  a 
linear  function  is  convex,  Lemma  3.5  yields 

W<P6^Pa<sp))  -  sA^3^(P3^P^(p))  <  s  Mi(X„^(Po\)). 

p  «  C(l)  and  s  >  0, 

and  similarly  for  P^(p) .  Then,  since  a  sum  of  convex  functions  is  convex  and 
since  the  square  of  a  nonnegative  convex  function  is  convex, 

(3.10)  W(Po^Pl(t.sp))  ^  s  g0a4x<Q0>-  P  ‘  C(l).  s.  t  ^  0. 

where 

<311)  '’Oi..x<«o)  -  <‘1<I1b«<’’o\)^  * 

Ve  factor  G(p)  as 


(3.12) 


G(p)  -  GqGj^(P) 
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where  Gq  Is  Independent  of  p  and  Gj^(p)  Is  linear  in  p.  Since  [L‘^P^(P)Gq]  and 
[Gj^(p)L*^)  are  linear  functions  of  p, 

convex  functions  of  p,  and  similarly  for  Pjj(p).  Therefore,  Lemma  3.5  and 
elementary  properties  of  yield 


(3.13)  ^  2  Pi(a„^(L*^Pl(t.p)Gosinwt) )  .  PiCa^CGiL*"^) ) . 


(Recall  Mi(*)  from  Definition  3.4.).  Proa 


p  <  C(l). 


(3.14)  2  (PgCp)  coswt  +  Pjj(P)  sinwt)  sin<i>t  - 

Pjj(p)  +  P^Cp)  sin2wt  •  Pjj(p)  cos2wt, 


it  follows  that 


(3.15)  2  <,„„<L-lpi(t,p)OosIiwt)  s  o„„(L-lpi,(p)Oo) 


Hence  (3.13),  (3.15)  and  Lemma  3.5  yield 

(316)  <  s2  ai„^(Qo)  «^2max<%>  •  P  ‘  ^(1).  s.  t  >  0. 

where 

(3-17)  ‘^laax<Qo)  "  ^l(‘^aax(L‘\Go)) 

and 

(3.18)  ‘'2max(Q0>  "  >*l(<^aax<(5ll'’’^^ > • 


Now  we  define 

(3.19)  si(Qo)  -  l/aax(«roBax(Qo).  (<^ia„(Qo)  •  ‘^2aax(Qo>)^^^»  • 


For  p  in  Che  interior  of  C(sj^(()q)),  it  follows  froa  (3. 10) -(3. 11)  chat 
Condition  3.1  holds  and  Ic  follows  froa  (3. 16) '(3. 18)  chat  Condition  (3.2) 
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holds.  Therefore,  we  have  the  following  theorem,  which  is  the  main  result  of 
the  paper. 

Theorem  3.6.  For  each  p  in  the  interior  of  C(sj^(Qq)),  P(t,p)  is  a  strict 
Lyapunov  matrix  for  A(t,p). 

From  (3. 3) -(3. 5),  it  follows  that  Pj^(t,p)  is  proportional  to  1/w  for 
large  w.  From  (3.11)  and  (3.17)  then,  it  follows  that  ‘^Omax^^^  ‘^Imax^^O^ 
are  asymptotically  proportional  to  1/w.  Since  q^„nv(Qn)  independent  of  w, 

<‘^lmax^%^  *  dominates  qQmax^Qo^  large  w,  so  that  Sj^(Qo)  is 

proportional  to  w^^^  for  large  w. 


4.  Numerical  Solution  of  the  Lyanpunov  Equations 

Eliminating  Pj3(p)  from  (3.4)  and  (3.5)  yields 

(4.1)  w2  P^(p)  +  +  Pa(p)A^  +  l^Jp^(p)AQ  - 

(C^(P)  Pq  +  Pq  G(P)). 

The  Bartels-Stewart  algorithm  [BSl]  for  solving  standard  Lyapunov  algebraic 
equations  can  be  generalized  in  the  following  way  to  solve  (4.1)  for  Pa(p). 

Let  U  be  a  real  unitary  matrix  such  that  U^AqU  -  Ag  where  Ag  has  quasi  Schur 
form.  Then  premultiplying  (4.1)  by  U^,  postaul tip lying  by  U  and  inserting  UU^ 
where  needed  yields 

(4.2)  J  Pg(p)  +  Aj^Pg(p)  +  Pg(p)Aj  +  2AjPg(p)Ag  - 

«  U'^IC^(P)  Pq  +  Pq  G(p)]U 

where  Pg(p)  -  U^Pg(p)U.  The  various  1x1,  1x2,  2x1  and  2x2  blocks  of  (4.2)  can 


be  solved  recursively  as  in  [BSl], 


8 


5 .  Example 

The  matrices  Aq  and  G(p)  in  (2.2)  are  the  following  4x4  matrices: 

(5.1)  Aq  -  I  0  I  I  G(p)  -  I  0  0  I 

l-Ko  -D  I  l-K^  0  I 

where 

(5.2)  Kq  -  |1  0|  D  -  .05  K<j  -  IPl  ^2^  ■ 

|0  4|  |P2  P3I 

Figure  1  shows  Sj^(I)  as  a  function  of  u  for  G(p)  factored  as  in  (3.12)  with 

(5.3)  Gq  -  |0|  (4x2),  Gi(p)  -  l-Ki  0|  (2x4), 

III 

and  for  G(p)  not  factored  (i.e.,  Gq  -  I  and  G2^(p)  -  G(p)).  That  Sj^(I)  is 
asymptotically  proportional  to  as  predicted  in  Section  3,  is  clear  from 

Figure  1. 

Perhaps  more  interesting  are  the  local  minima  at  w  >•  1,  2,  3  and  4.  We 
recall  a  classical  result  for  the  undamped  Mathieu  equation  (see  [NMl]  or 
other  standard  references):  a  parametric  excitation  of  frequency  twice  that  of 
the  nominal  system  makes  the  solution  to  the  equation  unstable.  Thus  the 
local  minima  at  u  -  2  and  u  -  4  (twice  the  natural  frequencies  of  the  nominal 
system)  might  be  expected.  Furthermore,  results  in  [NMl,  Chapter  5]  for  a 
i8ultl-degree*of’freedom  Mathieu  equation  indicate  additional  instabilities 
produced  by  parametric  excitation  frequencies  equal  to  s\ims  and/or  differences 
of  natural  frequancies  of  the  nominal  system.  (For  the  one-dagree*of 'freedom 
damped  Mathieu  equation  obtained  by  taking  Aq  and  G(p)  to  be  2x2  matrices  of 
the  forms  in  (S.l)  and  (S.2)  with  Kq  -  1  and  -  Pj^,  we  have  obtained  an 
S2^(I)  plot  similar  to  Figure  1  but  with  only  the  local  minimum  at  w  -  2.) 
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Figure  5.1.  ConiroUcr  Order  Variable 

Payload  Mass  —  0.1  a  Flexible  Link  Mass 
Gravity  Included 


Figure  5.5.  Conirollcr  Order  Variable 

Payload  Mass  *  0.1  x  Flexible  Link  Mass 
Gravity  Not  Included 
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Figure  5.6.  Controller  Order  Variable 

Payload  Mass  -  0.7  x  Rexible  Link  Mass 
Gravity  Not  Included 


